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PREFACE 

As  a  contribution  to  The  RAND  Corporation's  studies  in  system 
analysis  and  synthesis,  a  scheme  is  delineaced  in  this  report  for 
approximating  prescribed  system  characteristics  by  sums  of  exponentials 
or  by  rational  functions.  The  analytical  results  and  computational 
aids  which  are  derived  are  applicable  to  a  broad  class  of  system 
optimizations,  such  as  those  found  in  radar  filter  design;  numerical 
examples  are  provided  to  illustrate  their  application  to  practical 
and  important  design  problems.  The  contents  of  this  Memorandjm  should 
be  of  interest  to  the  A< r  Force  Systems  Command,  as  well  as  to  others 
concerned  with  numerical  methods  of  syftem  design  and  signal  analysis. 


SUMMARY 


In  many  system  design  problems,  the  representation  of  certain 
system  attrioutes  uust  be  in  terms  of  exponentially  damped  sinusoids 
or  of  rational  functions  in  order  to  be  physically  meaningful.  In 
this  Memorandum,  two  sets  of  orthonomal  elements  are  derived  which 
should  be  useful  for  such  approximation  problems.  One  set  constitutes 
a  basis  for  exponential  approximation  and  the  other  a  basis  for  ratiunal 
function  approximation. 

The  closure  properties  of  the  two-parameter  exponential  and 
rational  functions  are  examined  first.  Expressions  are  then  presented 
for  efficiently  determining  the  orthon  rmal  elements  of  each  basis. 
Important  characteristics  of  the  sets  are  also  discussed,  and  .pecial 
relations  among  the  coefficients  generating  the  bases  are  decuced. 

Once  the  general  relations  for  the  orthononnal  elements  are 
developed,  the  two  seta  are  applied  to  typical  approximation  problems 
encountered  in  iignal  processing  and  system  design.  The  computations 
involved  in  the  solut -on  of  these  problems  are  illustiated  in  the 
final  portion  of  the  study.  The  computer  programs  used  for  the  sample 
problems  are  described  in  the  appendices  and  should  be  helpful  in 


similar  design  situations. 
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I .  INTRODUCTION 

To  ensure  physical  realizability  in  the  synthesis  of  a  system,  it 
is  frequently  necessary  to  represent  certain  of  the  system's  r*: tributes 
by  approximants  other  than  algebraic  polynomials.  Such  is  the  case, 
for  instance,  in  the  determination  of  optimum  filters  for  smoothing 
and  prediction  where  it  is  convenient  to  approximate  empirical  or 
analytical  spectra  of  random  processes  by  readily  factorized  rational 
functions/  The  utility  of  exponential  functions  is  well  known  in  the 
design  of  linear  networks  for  a  prescribed  transient  response.  Simi¬ 
larly,  sums  of  exponentials  often  afford  the  most  suitable  approxima¬ 
tions  to  cross-correlation  measurements,  to  radioactive  decay  and  gas 
absorption  data,  to  mass  spectrographs,  and  to  ultracentrifuge  analysis 
curves . 

In  the  present  treatment  of  rational  and  exponential  function 
representations,  two  orthonormal  bases  for  exponential  and  rational 
function  approximations  are  derived.  The  bases  consist  of  two- parameter 
elements  which  provide  more  efficient  minimum  mean-square-error  approxi¬ 
mation  than  two  corresponding  one-parameter  sets  investigated  previously. 
After  the  closure  properties  of  the  two  orthonormal  bases  are  examined 
in  detail,  new  expressions  are  developed  tor  efficiently  generating  the 
individual  orthonormal  elements.  Several  relations  are  then  deduced 
which  connect  important  properties  of  each  basis.  Useful  Identities 
and  numerical  techniques  Involving  the  basis  coefficients  are  derived 

t 

The  term  "rational  function"  denotes  a  ratio  of  algebraic  poly¬ 
nomials  in  which  the  denominator  polynomial  is  not  generally  a  constant. 


(1) 


which  obviate  storage  of  either  the  fore  or  selected  values  of  the 
orthonoraal  approx iaants.  Finally,  several  numerical  examples  are 
provided  to  illustrate  both  algorithms  for  exponential  and  rationa 
function  approximation. 


1 
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II.  SCOPE  OF  INVESTIGATION 


In  2  previous  investigation,  the  linearly  independent  functions 


*„(t) 


-not 

e  0  <  t  <  * 


t  s  0 


(1) 


n  =  1,2, _ ;  Im(o)  =  0;  Re(o)  '< 


ind 


V  (tc)  =  - — 

fn  no  +  joj 


-ec  <  tc  < 


(2) 


n  =  1,2,...;  Im(o)  =  0;  Re(o)  >0 


were  examined  in  detail*  The  more  general  sets  comprised  of 


<p„(0 


0  <  t  .. 


t  i  0 


(3) 


n  «  1,2,...;  Im(o)  =  Im(P)  =  0;  Re(a)  >  0 


and 


*n^  no  +  j(u)-np) 


(4) 


n  =  1,2,...;  Im(o)  *  Im(|3)  =  0;  Re(o)  >  0 


will  be  treated  in  this  Memorandum.  The  primary  motivation  for  such 
a  study  is  chat  for  a  fixed  number  of  orthonormal  elements,  an  expan¬ 
sion  of  a  prescribed  function  in  terms  of  m  (t)  or  i  (w)  provides  a 

n  *  n 

minimum  mean-square-error  less  than  or  equal  to  that  for  q>^(t)  or 
^(«u)  ,  recpectively . 


t 

See  Tables  1-7  of  Ref.  1. 
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In  the  remaining  sections,  the  results  obtained  for  the  one- 
parameter  sets,  Eqs.  (1)  and  (2),  are  extended  to  {q>n(t)}  and 
and  numerical  schemes  ars  derived  for  simplifying  least- square  repre- 
sentatione  involving  {q>n(t)}  and  In  so  doing,  the  following 

topics  are  considered; 

o  Closure  of  and 

o  Determination  of  the  orthonormalized  setst 
m 

X<t)  =  T  Ymn  <P„(0,  Y  (t)  =  ~  [X  (t)  +  X*(t)] 
m  4^  mn  Tn  tn  /  £  m  in 

n«=  1 

ard 

m 

V"”  ■  I  '  „„  ♦„<“>  •  V")  "  75  &.<•>  +  ] 

n=  1 

o  Properties  of  X  ,  Y  ,  U  ,  V  ,  y  ,  and  \ 

r  m’  m  or  nr  Tmn’  mn 

o  Selection  of  a  and  0 
o  Computational  aspects  of  approximants 

M  M 

g(t)  «  Y  a  X  (t)  and  h(a>)  mV  b  U  (u>) 

'  4,  m  m  jLmm 

m=  1  m*=  1 

- 

The  symbol  *  denotes  complex  conjugate. 


III.  CLOSURE 


In  order  to  approximate  prescribed  functions  arbitrarily  closely 

by  linear  combinations  of  the  elements  cp  (t)  and  *  (<u)  ,  it  is  necessary 

to  verify  the  closure  of  the  sets  {cp  (t)}  and  {^(ui)}.  Since  the 

functions  <p  (t)  and  4  (uu)  are  related  through  a  linear  operation,  the 
n  Tn 

Fourier  transform,  it  will  be  possible  to  demonstrate  closure  of 

2  2 
{^  (u>)3  in  the  space  L  (-co,co)  once  closure  of  {cpn(t)j  in  L  (0,oo)  is 

shown. 

An  assemblage  of  functions  {vn(0}  of  integrable  square*  is  closed 
over  (a,b)  if  the  integral 

y(t)  cp*(t)  dt  =  0  n  =  1,2,...  (5) 

2 

implies  that  y(t)  C  L  (a,b)  vanishes  everywhere  in  (a,b)  except  perhaps 
over  a  set  of  measure  zero.  The  closure  property  derives  its  signifi¬ 
cance  from  a  theorem  that  states  that  a  set  of  functions  is  closed  if 

(2) 

and  only  if  it  is  complete.  '  Moreover,  if  a  set  {cpn(t)j  is  complete . 

2 

then  for  any  function  y(t)  c  L  (a,b)  and  any  positive  s,  there  is  a 
sum  function 

N 

5<t)  *  l  \  (6) 

n=  1 

such  that 

fb  2 

I  |y(t)  -  y(t)  p  dt  <  f  (7) 

w 

a 

r  2 

This  property  is  symbolized  by  cp^(t)  c  L  (a,b). 
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Equations  (5)  and  (6)  represent  a  special  case  of  the  general 

situation  in  which  a  closed  finite  or  infinite  system  of  elements 

{y^Cx)}  in  a  Hilbert  space  Y  permits  arbitrarily  close  approximation 
2 

(in  the  L  norm)  to  every  element  y(x)  c  Y  by  a  finite  linear  combina¬ 
tion  of  the  y^(x) .  To  establish  this  property  for  the  elements  cpn(t) 
in  Eq,  (3),  Szasz's  theorem*  can  be  invoked; 


A  necessary  and  sufficient  condition  for  closure 

2 

L  (0,1)  of  the  functions  x  ,  Re(\  )  >  -k , 

n  ‘ ' 

is  divergence  of 


V  1  +  2Re^n3 

Ji 1  +  M* 


(8) 


From  the  previous  definition  of  closure,  if  {x  }  is  not  closed 
2  2 

in  L  (0,1),  a  function  y(x)  CL  (0,1)  exists  which  is  orthogonal  to 

xn 

every  element  x  ;  that  is 


P*  2 

o  <  J  |y(x)  |  dx  < 


(9) 


and 


r1  K 

f  y(*)  * 

•  rv 


dx  -  0 


1  2 

1  ' 


(10) 


2 

Conversely,  if  the  set  {x  J  is  closed,  a  function  y(x)  CL  (0,1), 

which  is  not  identically  zero,  does  not  exist  so  that  Eqs.  (9)  and 

(10)  are  satisfied.  Accordingly,  with  the  change  of  variable  x  =  e* 
X 

in  the  set  [x  n],  the  conditions  given  in  Eqs.  (9)  and  (10)  become 


t 

See  Ref.  2,  pp.  32-36. 
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0  < 


(9  0 


and 

p°  t  C<1+U 

j  y(e  )e  n  dt  *  0  n  =  1,2,... 


(10  ') 


t  t/2 

With  the  further  transformation  y(e  )  e  =  y(t)  ,  the  last  two 
relations  give 


0  < 


dt  <  oo 


(9") 


and 


e  dt  -  0 


n  =  1,2, 


(10*) 


^  _  2 

Thus,  closure  of  the  functions  x  in  the  space  y(x)  c;  L  (0,1) 
is  equivalent  to  existence  of  a  function  y(t)  satisfying  Eqs.  (9*) 

and  (10*).  This  in  turn  is  equivalent  to  the  closure  of  the  functions 

t(^+X  )  2  ~ t(%+X  )  2 

e  in  L  (-«#,0)  ,  or  e  in  L  (0,®). 

In  the  application  of  interest,  where  cp^(t)  =  e 

satisfies  the  equation 


Xn  *  (-k  +  ns)  -  jng 

(11) 

n  =  1,2,...;  Im(a)  =  lm<0)  =  0;  Re(a)  >  0 


Consequently, 

Re[Xn]  =  +  na  >  n  =  1,2,...;  Re(cr)  >  0,  Im(a)  =  0  (12) 


as  required  by  Szasz's  theorem.  Substituting  this  value  of  X  in 

n 

Eq.  (8)  yields 
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y — i ±  2(-y  +  «>)  .  y  -  <»> 

^  1  +  (-^  +  no)  +  (n@)  n  (O  +8  )  -  no  +  5/4 


According  to  Szasz's  theorem,  closure  of  {?pn(t)  =  e  n^°  ^*3  in 
_  2 

the  space  y(t)  c  L  (0,«»)  rests  on  divergence  of  this  sum. 

Application  of  the  integral  test^  to  the  sum  in  Eq.  (13) 
indicates  that  divergence  of  the  infinite  series  depends  on  divergence 
of  the  integral 

f  _ 2q£  dL_ -  (14) 

Ji  i  +  (-Jj+co2  +  (es)2 


The  resultant  integration  is 


2c 


(oVw? 


7 TT  tan 


-l 


4a  +5p 


+  53 


(15a) 


2  2 
(o  +8  ) 


log  jV+p2)  §2  -  a§  +  5/4] 


or 


2a 


2  2  2  2  It 

(a+0Z)(4a+5pz)* 


tt  -1 

2  -  tan 


2  2 

2(o  +8  )-o 


2  2  k 

(40  +59  )* 


>- 


2  2 


2  2 

(a +8Z) 


log  [a  +8  -Of 5/4] 


(15b) 


+  lim  — ~T-  l°g  [(02+B2)  ?2  -  +  5/4l  “*  00 


2  2 
5-®  (o  +0  ) 


0  o 

Since  0  and  8  are  positive  and  bounded,  the  last  term  of  Eq.  (15b) 
is  unbounded.  Consequently,  the  series  of  Eq.  (13)  is  divergent,  and 
the  closure  conditions  of  Szasz's  theorem  are  satisfied  by  the  set 


-9 


f°ll°W8>  therefore,  that  {cpn(t)  =  e  n^a  is  also 

2 

complete  in  the  space  of  functions  y(t)  CL  (0,®). 

2 

The  closure  of  [^(ou)]  in  the  transform  space  D(uu)  CL  (-®,®) 

is  deducible  directly  from  a  lemma  of  Wiener's  on  invariance  of 

closure.’  His  lemma  states  that  the  closure  of  a  set  of  functions  is 

2 

preserved  in  any  linear  transformation  which  carries  the  whole  of  L 
into  itself  and  which  conserves  the  integral  of  the  squared  modulus  of 
each  function.  Quantitatively,  the  lemma  states  that: 


Given  a  linear  transformation  such  that  to  every  function 

2  2 
f(x)  CL  (a,b),  there  corresponds  a  g(y)  CL  (c,d),  if 

-  gjOO  and 

i.  c  f^(x)  -  c  gj (y) 

ii. "  f^x)  +  f 2 (x)  -  gx(y)  +  g2(y) 

iii.  |f  (x)|2  dx  =  f  |g  (y)  |2  dy  j  =  1,2,... 

a  J  c  J 

then  the  closure  properties  of  a  sequence  {f  (x)J  are  the  same 
as  those  of  (gj(y)}. 

The  transformation  which  carries  the  set  fm  (t)}  into  the  set 

Tn 

f*  (<»)}  is  the  Fourier  transform,  since 
n 


3  (<pn(t)}  ^  f  «pn(t)  e-jU)t  dt  =  f  e’lKHlv-Mn  dt 
-®  "  0 

_e-(nw-ju)-jn0)t  “  L 

no+jw-jng  ”  na+j(ui-n0)  “  ^n^ 

0 


See  Ref.  2,  pp.  28-30. 

1 1 

^The  symbol  -•  denotes  correspondence. 
3  denotes  the  Fourier  transform. 
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Because  the  Fourier  transform  is  a  linear  operator,  properties  (i) 

and  (ii)  of  Wiener's  lemma  are  satisfied.  Moreover,  Parseval's 
(4) 

theorem  indicates  that 

f  k(t>  I2  “  -  i  f  *•-  f  It„(2n*)i2  df  (17) 

"0  -JO  “ 00 

where  ou  =  2rrf.  Hence,  all  the  conditions  of  the  closure  invariance 

2 

lemma  are  fulfilled,  and  the  set  {i|i  (2nf)}  is  also  closed  in  L  (-<*>,«)  . 


-11- 


IV.  ORTHONORMALIZATION  OF  {cp  (t)J  and  {^(uu)} 


The  closure  properties  discussed  in  Section  III  ensure  completeness 

2 

of  ftp  (t)j  in  the  space  of  squared- integrable  functions  d(t)  cL  (0,») 

2 

and  completeness  of  Un(“>)  J  in  D(ou)  c  L  (  —  co j oo)  ,  Consequently,  the 
theory  of  Hilbert  spaces  can  be  applied  to  two  further  items  of  com¬ 
putational  importance:  orthonormalization  of  {tp^Ct)}  and 

2 

and  representations  of  functions  d(t)  CL  (0,»)  by  sums  of  elements 
cpn(t)  and  tn(u>)  . 

The  task  of  orthonormalizing  the  set  {cpn(0]  (or,  analogously, 

the  set  f  \jrn (uo)  J )  entails  finding  coefficients  y  in  the  functions 

X(t), 

m 


x  (t)  =  y  y  ®  (t) 

m  '  Tmn  rn 


m  =  1,2,... 


such  that 


(*®  * 

J  Xr(t)  Xfl(t)  dt  =  6i 


7.he  existence  of  these  constants  y  is  guaranteed  by  Theorem  IV.  1 

Tmn 

of  Ref.  5: 


Let  cp-i  > *$2 » •  •  •  *  a  finite  or  infinite  sequence  of  elements 
such  that  any  finite  number  of  elements  tp^  • -<PK  4sre 
linearly  Independent.  Then  constants 


V21  Y22 


Y3 1  y32  y33 


The  symbol  6  is  the  Kronecki r  delta:  6  -  1  for  r  »  a:  0  -  0, 

.  .  ra  rs  rs 

otherwise . 
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can  be  found  auch  that  the  elements 

xi  ■  ''n  n 

X2  ‘  *21  *1  +  Y22  *2 

xj  -  Yn  *!  +  Y32  92  +  V33  9j 


are  orthonormal. 

The  proof  of  this  theorem  provides  an  iterative  algorithm,  the  Gram- 

Schmidt  orthogonalization  procedure,  for  actually  determining  the 

functions  X  (t).  The  recursion  is  given  by  the  following: 
m 

”  «PX  and  \  *  x1/llac1H 

X2  *  *2  *  r'*2>Xi>  Xi  and  X2  m  x2/Hx2JI 

!  (20)’ 

m 

X»1  ■  Vi  •  I  <Vi-V  \  *,’d  X»1  •  Vi^Vi11 

k-1 


Although  the  functions  X  (t)  can  be  determined  by  this  iterative 

tn 

procedure,  Eqs,  (20)  provide  only  an  implicit  expression  for  the 
coefficients  v  generating  fX  (t)l.  An  additional  drawback  to  the 
above  scheme  is  that  the  recursive  evaluation  of  the  basis  elements 
is  exceedingly  laborious  .M 


*The  notation  jjf||  denotes  jj  f(t) 
•  L*0 
f  f(t)  8*(t)  dt.  Thus,  ||f||2  -  (f,f>. 

*  A 


t  I 


See  Ref.  1,  pp.  12-14. 


and  (f,g)  symbolizes 
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To  circumvent  these  difficulties,  it  is  necessary  to  abandon 

the  Gram-Schmidt  approach  and  to  consider  obtaining  the  ymn  from 

methods  based  on  the  set  {^(<d)}  derived  from  the  Fourier  transform 

of  cpn(t).  With  the  aid  of  Parseval's  theorem  and  other  key  results 

from  the  theory  of  functions  of  a  complex  variable,  it  will  be  possible 

to  achieve  an  explicit  relation  for  the  coefficients  v  >  as  well  as 

'mn 

for  the  transform  parameters  X  in  the  equation 

mn 


where 


U  W>  -  Y  X  4i  (u>) 
m  mn  Tn 


OD 

J  ur(«j)  U*(ou)  duj  =  6j 


Application  of  the  Fourier  transform,  Eq,  (16),  to  the  elements 
cp^(t)  and  i^(<d)  of  Eqs,  (3)  and  (4)  reveals  that 


J  U  (u>)  e^1  du>  =  4l  J  Y  X  1  (is)  e^  du> 

2tt  J  ru  2n  J  A  mn  Tn 

-00  -  00  , 

n=l 

m  m 

-  y  jTp  f  t  (uj)  e^  duu  =  y  X  cp  (t)  =  X  (t) 

A  J  Tn  mn  Tn  m 

■  OD 


where  m  is  finite.  Hence,  X  (t)  and  U  (u>)  are  Fourier  transform  pairs. 

m  m 

Since  U  (id)  is  assumed  to  be  orthonormal,  Parseval's  theorem  leads  to 
m 

*  p®  *  '••ff 

U  (id)  U  (id)  du)  =  2tt  X  (t)  X  (t)  dt  =  6  (24) 

w  m  n  J  m  n  mn 


Consequently,  if  the  coefficients  v  in  X  (t)  are  chosen  as 

mn  m 


'i 


then  the  set  {X  (t)}  will  also  be  orthcnoraal ,  for 

9 


C  *  2tt 
cm 


I  I  *.«  V0  [  I  4  4"  1 

■*J-1  J  [fc-i  J 


■  1  ^  x=»  !  I  ^  4  "v(t) 

*°  n*  1  |^k=*  1 

’  to  1  f  n 

00  _ 

'  J  1  V.  *n(t)  I  4  4°  dt 

J  jk=l 

oo 

•=  r  x  (t>  x*(o  dt 

J  m  n 


Once  the  orthonoraal  elements  (U^C®)}  are  found,  therefore,  the 
orrhonorvnai  set  [X  (t)]  is  aiso  uniquely  determined. 

9 

The  condition  of  Grthononualitv  and  the  structure  of  f  (x) 

■  n 

suggest  a  scheme  for  finding  the  parameters  \  of  U  («) .  The 

mn  m 

procedu  e  is  based  on  two  results  of  analytic  function  theory,  and 
is  an  extension  of  the  application  of  the  theorems  found  in  Refs.  1 
and  5. 

By  observing  the  pole  patterns  of  the  rationalized  functions 

U  (®) ,  m  =  1,2,...,  it  is  possible  to  compute  the  coefficients  \ 
m  mn 

which  determine  the  m^  orthonormal  function  U  («)  .  In  U.  (®)  ,  for 

m  I 


example , 


U1W  *  „  + 


?  pole  is  located  at  ji  =  -o  -  jg.  In  L'^Ob),  where 


l'2(«) 


CH- j  (sc-9)  2w-j  (0.-26) 


-(2a.2i+^22^  +  2 l-a;^'21+i'22^  ~  ^^21+*z2^ 

;«-j(x-g)j  C2w-j(x-2g)2 


poles  are  located  at  jx  =  -w-jp  and  jx  =  -2ofj2g.  In  general,  for 


U  (x)  =  - — - +  - — - +  •  -  •  +  - — - 

a  CH-j(x-g)  2of  j  (x-2g)  ao*- j  (x-ag) 


equally  spaced  left-half  s-plane  (LHP)  poles  are  located  at  ju  =  -0+jg, 

* 

-2:»-j2g,  ....  -acH-jag.  Siailarlv,  the  conjugate  function  U  (u?)  has 

E 

equally  spaced  right-half  s-piane  (RHP)  poles  at  jx  =  c+jg,  2c+2jg, 
ao+jap.  Therefore,  if  U  (x)  is  expressed  as 


V'i(x)  =  — .  ^  ,  a i  =  r.oraalizing  constant  (30) 


then  to  orthogonalize  U^(x)  and  U^(x) ,  i.e.,  to  ensure  that 


j  i'2(x)  l^x)  ix  = 


U2(x)  aust  be  chosen  as 

a2  [ff-j(«f-3)j 

I'  (x)  =  — —  ~  FT ~  ,  a  =  normalizing  constant  (32) 

2'  Lcw-j(x-p)j  i2c*-j(x-2g)]  ’2 


In  this  way,  the  zero  of  U  (*)  cancels  the  RHP  pole  of  U  (to)  in  the 

^  1 

* 

product  l’^,  and  the  integrand  of  Eq.  (31)  becomes  analytic  in  the 
*  (6) 

RHP.  Since  U  U  also  satisfies  Jordan's  lemma,  the  Cauchy  integral 


theorem''  can  be  applied  to  the  integral  in  Eq.  (31)  (with  the  contour 
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of  integration  taken  in  the  clockwise  sense  as  the  jco-axis  and  infinite 


RHP  semicircle)  to  establish  that 


J  U2(u)>  U^(tu)  duu  = 


* 

alfl2 


ds 


[afa-jp]  [w-2a-j2g] 


=  0 


(33) 


Similarly,  if  U^(uj)  is  selected  as 

a3[o-j(u)-e)]  [2o-j(u)-2p)] 

U3^  "  [w-j((U-0)]  [2of  j  (u>-20)  ]  [3w-j(tu-3e)] 


a^  =  normalizing  (34) 
constant 


then  both  U2  and  are  analytic  in  the  RHP,  and  consequently 


both  equalities 


* 

r®  *  -i  J*  a2a3  d8 

J  °3«  V«»  d“  '  ~f.  [W2a.i28]  *  ° 

-  cd  -  J  » 


(35) 


and 


r”  *  .i  a*a3  [-s+2ofJ2p]  ds 

J  °3W  *• '  1  9  ,  [»»-isn*2a-2jj]  |>3o-J3B] 

-00  -  J® 


=  0 


(36) 


are  guaranteed  by  Jordan’s  lemma  and  the  Cauchy  integral  theorem. 
Extending  the  above  sequence  to  the  m^  element  of  the  orthogonal 
basis,  it  is  observed  that  the  general  form  of  U^iu)  must  be 


1 

<W-jOu-0) 


U  (uj)  =  i 
m 

a  [a-j(iu-g)]  [2a- j (u)-2g)  ]•  •  •  [(ra-l)a- j  (u>-mfH-g)  ] 
m _ _ _ _ 

[c+  '  'u)-g)  J  [2afj  (u>-20)  ]•  •  •  [maf  j  (w-mg)  ] 


(37) 


m 


2  3 

*"  >  -*  I  *  *  * 


where  a  is  a  normalizing  constant, 
m 

The  constants  a  are  readily  determined  by  requiring  that  the 
m 

set  [Um(u))J  be  orthonorraal;  i.e., 
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00 

P  * 

U  (<d)  U  (id)  do)  = 
J  m  m 


Substituting  Eq.  (37'i  for  U  (<d)  in  the  normality  condition  Eq.  (38) 

m 

gives 

r“  ,  %  Iam|  2[a"J  (tt>- 3)  ][2o-j  ((D-2g)  ]. . . [(m-l)o-j  ((D-mgf  p)  ] 

J.  "  j  [of  j  ((D-g)  ][2of j  (uj-23)  ]. . .  [mof  j  (<D-mg)  ] 


LCH-.1  (u)-B)  11  2ch-  j  ((P-2p)  1...I  (m-l)gfi((D-n 
[a- j  ((D-g)  ][2a- j  (id-28)  ]. . .  [ma-  j  (u>-mp)  ] 


00 

■  J  J 


a  |"  d(juj) 

_ m _ 

[mofj  (iD-mg)  ][(tKJ-  j  (iD-mg)  ] 


=  2nj 


j  mCTfjmp~j(D 


I  j(D=-nmjtnp 


=  1,  m  -  1,2,.., 


Since  the  U  (id)  are  already  orthogonal,  the  normality  condition  Eq.  (38) 
m 

can  be  satisfied  by  choosing  the  a  real  in  Eq.  (39)  and  by  assigning 

m 

to  a  the  value 
m 


m  =  1,2,, 


In  view  of  Eq.  (37),  the  orthonorrnal  basis  functions  can  be  finally 
written  as 


2  1 

TT  Oh  j  (lD“  p)  ’ 


m  =  1 


U  ((D) 
tn 


_  ft  [ no- J  ((D-ng)  ] 
mo  n=  l 


n  m 


m  =  2,3,, 


ff  [n<H  j  ((D-ng)  ] 
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The  coefficients  \  in  the  series  representation  of  U  (u>)  , 
mn  tn 

Eq.  (21),  can  now  be  related  to  Eq.  (41)  by 


/w  [q-j(ui-B)  1[2a-1  (m-26)  ]• 
V  n  [ofj  (a‘-p)  ][2ofj  (o>-2p)  ]• 


•  •[(m-l)CT-i(u)-m&fB)  1 

•  •  (iD-mg)  ] 


■  l 


mk 


kcH-j  (ou-kp) 


m  =  2,3,... 


k=l 


(42) 


Multiplying  both  sides  of  Eq.  (42)  by  nof  j  (u)-n0)  ,  1  £  n  £  m,  gives 


/no  [ngf  1  (w-nfl)  Ifo-Km-e)  ]F2ct-  j  (m-2B)  1-  •  •  f(m-  l)g- 1  (ou-m&fB)  ] 
V  n  [cH-j  (io-p)  ][2w-j  (ui-28)  ]•  •  •  [nw-j  (uj-np)  .  [mw-j  (u)-mp)  ] 


m 


k=l 


ng+  1(ui-nB) 
ko+j  (uJ-kp) 


m  =  2,3 , . . . ;  1  *  n  £  m 


(43) 


Since  Eq.  (43)  must  be  valid  for  all  tu,  it  must  also  be  an  identity 

for  u)  ■  npfjna,  or  jiu  =  -nofjn8.  Accordingly,  with  ju)  =  -nofjnp, 

Eq.  (43)  reduces  to 

m- 1 


fj  [c(rH-  r)  -  j8(n-  r)  ] 


J]  [-a(n-r)+jp(n-r)  ]  jj  [-o(n-r)+jp(n-r) ] 
r=  1  r=iw-l 


so  that 
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n  -  m  =  1 


mO  r=l 


If  [-Kf^)] 


If  [r-n] 


m  =  2,3,...;  l^n^m 


m  =  1,2,...;  n>m 


In  view  of  Eq.  (25)  ,  v  is  immediately  determined  for  the  ortho- 

Tmn 

normal  set  {X  (t)}  as  v  =  •'fzn  \ 

m  'mn  v  mn 

The  preceding  derivation  of  the  orthonormal  basis  (U^io)}  permits 

a  simple  determination  of  the  allied  real-valued  functions  of  uu 

m 

V»<“)  ■  J7  K™  +  Um(“)]  ’  /2  I  t\.„  +  C  ^  (46) 


Assurance  of  orthonormality  for  the  set  (V  (uu) }  follows  from  the 

m 


relations 


j  V  (is)  V*(is)  dcu  =  f  [U  +  U*J  [U*  +  U  ] 
J  m  n  2  J  1  m  mJ  n  nJ 


f  V  (ou)  V*  da>  =  f  U  U*  dis  +  -J  f  U*  U  dis  +  -J  f  U  U 
Jm  n  2  J  m  n  2  J  m  n  2Jmn 

-  00  -00  -00  -00 

< 

+  ”  f  U*  U*  dis 
*.  J  m  n 


For  the  functions  ^  (is)  defined  in  Eq.  (2),  the  result  for  \ 

n  _  mn 

with  0=0  agrees  with  the  quantity-/—  a  derived  in  Ref.  1. 

V  tt  mn 
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Since  the  U  are  orthononr.al  (Eqs.  (33-38))  the  first  two  integrals 
on  tne  right  side  of  Eq.  (48)  yield  the  Kronecker  delta 


f  V  V*  dm  *  6  +  f  U  U  dm  +  f  U*  U*  diu 

Jmn  mn  J  m  n  2  J  m  n 


(49) 


It  is  evident  from  Eq.  (37)  for  U  (uj)  that  the  integrand  U  U  in  the 

m  m  n 

first  integral  of  Eq.  (49)  contains  only  LHP  pole3.  Since  Jordan's 

lemma'  is  satisfied  for  U  U  ,  the  first  integral  of  Eq.  (49)  may  be 

m  n 

evaluated  over  the  infinite  RHP  semicircle  with  diameter  on  the  jui-axis. 

Since  Um  is  also  holomorphic  in  this  RHP  region,  Cauchy's  integral 

theorem**  guarantees  that  the  integral  of  U  D  vanishes.  Similarly, 

m  n 

*  it 

with  U  U  holomorphic  in  the  entire  LHP,  it  can  be  argued  that  the 
m  n 

second  integral  of  Eq.  (49)  is  also  zero.  Consequently,  Eq.  (49) 
becomes 

CD 


I 


V  (uj)  V  (»)  do)  «  6 
m  n  mn 


(50) 


-CD 


The  rational  function  expression  of  V  may  be  written  with  the 


m 


aid  of  Eqs.  (46),  (21),  and  (4)  as 


ra 


,  _  (X  +X*  )  no  +  j(-X+\*  )(uJ-n&) 

„  ,  ,  1  -r  am  mn  mn  mn 

Vm(U,)  “  75  )  - 2 - 2 

72  tl  (W)  +  («-n0) 


(51) 


or 


v  (IB) 
m 


n»l 


3  *  '  t 

inn  •  v  4  g 

s  +  (tu-ng) 


(no)' 


n=  1 


t 

See  Ref.  6. 
l  * 

See  Ref.  7. 
ft* 

For  8  ■  0,  X  is  real  and  b  ■  0.  When  /2a  is  also  an 
H  *  mn  mn  mn 

integer,  m,n  =•  1,2,...,  V  (<i>)  takes  the  form  of  the  envelope-delay 
components  examined  in  Ref.  1. 
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Substituting  these  relations  into  the  expression  for  V  (m)  , 


’  A  [u«w  +  u>>]  *  A 


m  m 

V  X  *  (ou)  +  7  X*  if*  (uu) 
mn  Tn  mn  Tmn 

n=  1  n=  1 


(46) 


it  follows  that  the  Fourier  transform  of  V  (u>)  is 

m 


M"  V“> 


j«ut  .  i_ 
e  da,.- 


m 


m 


y  x  9  (t)  +  y  x  9  (-0 

mn  “n  mr.  Tn 


n=  1 


n=  1 


(58> 


Consequently,  Parsevsl's  theorem  and  the  orthonormality  relation, 


Eq.  (50),  for  V  (w)  allow  Eq.  (50)  to  be  written  as 


m 


®Y(t)Y(t)  oo  *  , 

I  7 TT  T^T"  dt  -  J.  Vm(u,)  Vn(U,)  di"  =  IS  6, 


(59) 


so  that 


r”  * 

Y  (t)  Y  (t)  dt  -  6 
J  m  '  n  7  mn 


(60) 


From  Eqa.  (55),  (53),  (45),  and  (3),  the  orthonormal  basis  functions 
Y  (t)  can  be  finally  expressed  as  the  complex-valued  functions  of  t 


m 


l 


X  e 
mn 


-n(o-jg)  t 


Ym(t)  =  < 
m 


n«  1 


m 


/TT  y  X  e 
L  mn 


*  +n(<*fjB)t 


t  >  0 


m  =  1,2,, 


t  <  0 


(61) 


n»  1 


or,  expressing  Xmn  in  polar  form, 


-23 


' /n  I 


A.  e 
1  mn  1 


-na|tj  I 


n=  1 


cos 


["fit  +  |7|  argCX^)] 


+  Ttf 


j  sin 


n$  t I  +  ar 


g(\  )~ 

mn 


f  m  =  1,2, , . .  V  t 


(62) 


The  key  relations  established  for  the  orthonormal  bases  X  (t) 

m 

and  U  (uu)  and  the  equations  determining  the  orthonormal  basis  coef- 
m 

ficients  \  are  summarized  in  Table  1. 
mn 


EQUIVALENT  REPRESENTATIONS  0 


5  8 


e  % 
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6  U-<1  n 
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G  L-'vi  l|. 


a  1^4’ 
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•H 

4) 

T3 

a 

CU 

o 

0) 

U 

y 

g 

CO 

o 

M-C 

c 

CO 

■*-C 

c 

CO 

’X3 

01 

■U 

> 

•r* 

u 

4) 

<u 

*H 

T3 

p 

3 

0) 

0 

u 

h 

CO 

0) 

c 

CO 

g 

n 

3 

Q) 

r— < 

> 

o 

c 

u 

c 

41 

e 

X 

r< 

4J 

4) 

CO 

x 

<u 

4J 

AJ 

* 

o 

c 

c 

•r4 

4) 

T3 

CO 

c  , 

* 

0 

1  .  ^ 

rt>  oo 

4J 

\£> 

CO 

•a 

d 

o; 

co  ■ 

u* 

-  u 

0 

u 

i 

O 

u 

4-1  X 

CO 

U 

U  T) 

00 

0)  4) 

CO 

a,  C 

1  r-A 

o 

The 

E  41 

C  "O 

0 

<0 

4-1  UJ 

CO  *rA 

c 

• 

<fl  c 

l 

W  g 

w  Q 

1 

••• 

-25- 


V.  SPECIAL  PR01LRTIES  OF  THE  BASIS  COEFFICIENTS  \  AND 
_ mn 

THE  ORTHONORMAL  ELEMENTS  X  (t)  ,  U  (uj) 
_ m  m 


The  preceding  analysis  culminates  in  Eq,  (45)  for  the  basis  coef¬ 
ficients  \  which  generate  the  orthonormal  elements  U  (u))  ,  V  (uo)  , 
mn  mm 

X  (t) ,  and  Y  (t).  In  order  to  detect  errors  in  the  computation  of 
m  m 

these  \  ,  it  is  de3irable  to  provide  a  check-sum  relation  for  the 

mn 

\  analogous  to  the  expression  for  the  coefficients  a  of  Ref.  1. 
mn  mn 

It  is  shown  in  Ref.  1  that  for  0  =  0  an  identity 

m 

V  &  =  (-1)"*1  m  =  1,2,...  (63)’ 

L>  mn 

n=l 

exists  in  which  the  a  ate  given  by 

mn 


See  Ref.  1,  p.  22. 

1 1 

See  Ref.  1,  p.  21. 

1 1 » 

See  Ref.  1,  Table  1. 
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For  the  »  derived  in  this  Memorandum .  it  is  similar iv  demonstrable 
an 


V  w?  I  Nan  *  1  a=”  (_1) 


a  =  1.2, 


n=  1  n=  1 


where  a  is  defined  as 
ten 


arin  'U  ac  'an  i  a  • 


ft  [-«] 


!  J;  Cr'": 


P.  =  E  =  1 


f“  s  2  3  *  1  1  r  i  » 

“*  4*|Jf  M  *•  **  *~ 


=  =  1,2,  —  ;  n  >  = 


For  the  special  case  in  which  5=0,  it  is  clear  that  =  a  , 

En  sn 

and  the  validity  of  Eq.  (67)  follows  directly  from  the  identity 
Eq.  (63).*  *  In  general,  for  arbitrary  real  g,  it  must  be  established 


or  tnat 


n=  1  n- 1 

E- 1 

«  jj  (n  +  rz) 

Z  «,/ 

n-1  J]  ^  "  r*^ 


^  fth«j 


-  (-d^1 


«  c-D 


l  ■ —  / 

The  prime  used  in  the  product  notation^  (r-n)  sign! 


ties  that 


r  ranges  only  over  the  integers  1,2, _ ,n-l  ,:-h-1  , . . .  ,»  so  that  the 

result  is  never  zero, 
t » 

See  Ref.  1,  p.  20,  Eq.  (64). 


••r.e.e,  for  cenvf  nier.ce .  ;  f  coeplex  variable  r  is  defined  as 


z  = 


=  -  rg 

(£) 

Free.  me  de::r::;c?.  or  ::.e  Stirhr.g  refers  of  :he  first  find, 
z  possible  to  vrite  the  factorial  function  relation 


(71) 


(x-I)(x-2)  ---  >  -  (:«-!) I  s  ; 


\  «a=j  _=*  *■ 


O'  X 


(72) 


(x-:)  (x-i) 


r“ 

•  )  s-  x~  ‘  c-ir 


(23) 


Hence .  the  numerators  of  Eq.  (~i;  car.  be  vritter.  as 

_  t 

fj  (r-  -  rr) 


V  /n  \  <-(*)  AiN  *  ,  . 

n  ij  v =  ^  >  sn  i*y  (‘° 


ana  the  left  sice  of  Eq.  (7C-)  becooes 


T~ 

L 


fl,  <— '  Z 

ATT 


(r-r.)  n»i  jj  (r-r.)  k= 


1  s-  ‘j;  (“'  =  1  *r 


(75) 


Ir.  addition,  the  denominators  o:  Ec.  ( 7 C)  can  be  expressed  as 


ft  (r-r.)  =  (-1  r~-J?  (r-r)  =  (I) 


ri(n-n): 


:  S  r  >  n 


(76) 


Jf  (r-r.)  =  (-1) 


r^i  nljn^l  _ 


r=  i 


(-D 


f~  \ 

nW 


(77) 


sc  that  Eqs.  (69),  (7C)  ,  and  (75)  can  be  fomulated  as 
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l 


ilii 


«rt- 1 


ta  m 


on  b: 


l  l  <-»*"  ©  «-k 


n- ; 


txl  k-0 


From  this  last  relation,  it  is  clear  that  Eq.  (69)  can  be 
established  as  an  identity  if  it  can  be  verified  that 


m  o 


£  )  <-!)W"  Q  s<k>  „k  z“'k  .  (-l)Wl 


k*0  n=l 


or,  equivalently,  that 


“  <-i)k  s(k) 

l—**- 

k*0 


m 


I  <-»*  C) 


n=  1 


a-k 


-1  =  0  ta  =  1,2,. 


Since  Eq.  (80)  is  a  polynomial  of  degre>  a  in  z,  it  follows  front  the 

.(9) 


fundamental  theoreo  of  algebra^  7  that  Eq.  (80)  has  only  a  values  of 
r  which  make  the  left  side  equal  to  zero.  In  order  for  Eq.  (80)  to 
be  generally  valid,  therefore,  each  coefficient  of  zr,  r  =  1,2, _ ,m 


must  be  zero  and  the  coefficient  of  z^  trust  be  unity  (since  S^^/tn'. 

*  tn 


is  also  nonzero  for  k  *  1,2, _ ,m;  tn  =  1,2,...).  Consequently,  the 

proof  of  Eq.  (69)  rests  on  demonstrating  that 


<-l)-k  s("~k)  * 


l  <-»"  0 ■ 


a-k 


n=l 


0 


k  =  1,2,, 


k  =  0 


or,  since  S 


(0) 


o  s(E) 


t00 


1,  and  Sv  7  ^  0,  k  -  1,2,..., a,  tn  =  1,2,. 
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l  <-»"  (:)  - 1  <-»■  o 


k  =  1,2,. ...m-l 


(-1)  m!  k  =  0 


In  order  to  verify  Eq.  (82),  the  following  identity  is  utilized: 


Amf(x)  =  Y  ^_1)m  U  (”)  f(x+n) 


where  A  denotes  the  forward  difference  operator.  With  the  choice 


f(x)  =  x 


k  =  0,1,. . . ,m;  m  =  1,2,... 


Eq.  (83)  becomes 


.m  m-k 
A  x 


-  I  0  <*■>' 


It  is  possible,  now,  to  conclude  the  demonstration  by  applying  to 
Eq.  (85)  the  fundamental  theorem  of  difference  calculus: 

The  n1”*1  difference  of  a  polynomial  of  degree  n 


yW  -  y 


a .  xJ 


a  -i  0 
n 


th 

is  a  constant,  a  n1.  h,  and  the  (m-l)  '  difference  is 

n 

equal  to  zero.  The  first  forward  difference  is  defined 

as  Ay(x)  s  y(x+h)  -v(x) ;  the  second  forward  difference, 
2 

as  A  y(x)  «  Ay(x+h)  -  Ay(x);  etc.,  and  h  is  a  constant. 


Accordingly,  for  k  =  l,2,...,ra-l 


m  m-k 
Ax  =0 


k  =  1,2,... ,m- 1 


See  Ref.  8,  p.  823. 
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so  that  Eq.  (85)  becomes 


l  0  (W">' 


k  =  1,2,... ,m- 1 


For  the  particular  choice  x  =  0,  Eq.  (87)  becomes 


l  <-‘>r  0  -“'k  ■ 0 


k=  1 ,2 , . . . ,m- 1 


which  verifies  Eq.  (82)  or  Eq.  (81)  for  all  k  =  l,2,...,m  except 
k  «  0 .  With  k  =  0,  h  =  1,  Eq.  (85)  may  be  simplified  through  the 
finite  difference  theorem  to 


.mm  ,  V  /  i\m~n  /m\  /  , 

A  x  =  ml  =  )  (-1)  (  )  (x+n) 


or,  with  x  =  0 


1  (-1)"  (")  n"  -  (-1)"  «! 


Since  this  last  relation  is  recognized  as  Eq,  (82)  with  k  =  0,  Eq.  (82) 
and,  therefore,  Eq.  (69)  are  established  as  identities  for  m  =  1,2,...; 
k  =  0,1,..., m.  From  Eqs.  (68)  and  (69),  it  is  clear  that  the  coef¬ 
ficients  \  satisfy  the  check-sum  expression 
mn 


m 

l  V.n  ' 


m  =  1,2, . 


TVo  additional  derivations  simplify  the  calculation  of  the  initial 

values  of  the  basis  elements  X  (t)  and  U  (m) ,  as  well  as  the  evaluation 

m  m 

of  the  integrals  of  these  functions.  In  order  to  get  these  relations, 
it  is  necessary  to  verify  the  identity 


-31- 
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Equation  (98)  follows  immediately  from  the  finite  difference 
relation  given  in  Eq.  (83): 

m 

A™  f(x)  -  £  (-l)”"11  (”)  f(x+n) 

n«0 

If  f(x)  is  chosen  as 

f(x)  =  xk_1  (99) 

then  Eq.  (83)  states  that 

m 

Am  k-1  V1  /  *\®~n  /  vk-1  ,,nn\ 

Ax  =  £  (-1)  [nJ  (x+n)  (100) 

n=0 

The  previously  cited  theorem  of  difference  calculus  provides  the 
relation 

A”  xk_1  =0  k  =  2 ,3 , . . . ,m  (101) 

Consequently,  Eq.  (100)  becomes 
m 

V  ,  ,.ra-n  /m\  ,  .k-1  „  ,  „  „  . 

)  \n)  +  =  0  k  =  2»3 . m  (102) 

n=0 


With  the  choice  x  =  0,  this  simplifies  to 


m 


n=0 


or,  since  the  first  term  of  the  summand  is  zero 


m 


n*=  1 


,m 


This  last  equation  is  recognized  as  Eq.  (98) . 

It  remains  to  demonstrate  the  validity  of  Eq.  (97).  Since 
sW  is  given  by 


(103) 


(104) 
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S(1)  =  (- l)m"  1  (m-  1)  ! 
m 


m  =  1,2, 


(105) 


Eq.  (97)  becomes 


m 


I  <-»■  © 


(106) 


n- 1 

or,  adding  and  subtracting  the  term  for  n  =  0, 
m 

I  (‘1>n  (")  =  °  m 

n=0 

But  Eq.  (107)  follows  directly  from  the  binomial  expansion^^ 


(107) 


m 


o  ♦  x)" .  I  © 

n=0 

Thus,  with  x  =  -1,  Eq.  (108)  yields 
m 

°  -  I  ©  <-»" 


m  =  1,2, .. . 


(108) 


m  -  1,2,.  . 


(109) 


n=0 


Since  this  is  precisely  Eq.  (107),  Eq.  (97)  is  verified,  and  thereby, 
the  identity  Eq.  (92).  Finally,  utilizing  Eqs.  (92)  and  (68),  it 
follows  that 


y  -y- 


i 


(110) 


n=  1 


The  identities  given  by  Eqs.  (91)  and  (110)  can  be  used  to  evaluate 

X  (0)  ,  U  (0)  ,  f  X  (t)  dt,  and  f°°  U  (m)  du).  From  Eq9.  (18)  and  (25), 
m  m  mom 

U  -  oo 

X  (t)  becomes 
01 


See  Ref.  8,  p.  824. 


X  (t)  =  Y  Jl n  X  cp  (t) 

m  '  ^  mn  Tn 

n=  1 

m  =  1,2,... 

(Ill) 

This  can  also  be  expressed,  by  Eq.  (3)  for 

cpn(t)  ,  as 

m 

xco-vsy*  e-n(°-j8)t 

m  L  mn 

n=  1 

tn  =  1,2,...;  t  >  0 

(112) 

Hence,  X  (0)  becomes 
m 

m 

X  (0)  =  t/2tt  y  X  =  ('l)'nf  1 
m'  v  u  mn 

n=  1 

m  =  1,2,... 

(113) 

Similarly,  using  Eqs.  (4)  and  (21;,  Um(w) 

can  be  written  as 

m 

=  Y  mn  tna  +  .Kw-^)-) 

n=l 

m  =  1,2,...;  |w|  <  «> 

(U4) 

so  that 

,  ,  1  r  \»„ 

V°>  •  o  -  JB  I  ~ 

n=l 

m  =  1,2, . .  . 

(115) 

or,  in  view  of  Eq.  (110) 

r—  ni-  1 

U  (0)  -  Jc  <ot  -ua— 

m  VtnTT  (a  -  jg) 

m  =  1,2, .. . 

(116) 

The  areas  under  U  (tn)  and  X  (t)  can  be  ascertained  in  a  similar 
m  m 

fashion.  From  Eq.  (23) 

m 

f  U  («)  eju,t  dii)  =  X  (t)  =  y  X  cp  (0 
2tt  J  m  m  ^  mn  n 

"®  n=l 

If  cp  (t)  i-8  substituted  in  Eq.  (23)  ,  it  follows  that 
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00 

-T-  r  u  (cu)  dm  =  y  \  cp  (o)  =  y 

2n  J  tn  Lj  mn  Yn  L 


mn 


n=  1  n= 1 

Consequently,  replacing  the  sum  of  the  \  with  the  right  side  of 


mn 


Eq.  (91),  the  area  under  U  (ou)  becomes 


(117) 


OD 

j  U  (ou)  dou  =  (-l)™*"  Vanina 
J  m 


m  =  1,2, .. . 


(118) 


Since  Eq.  (23)  relates  U  (ou)  and  X  (t)  as  Fourier  transform  pairs, 

m  m 

it  is  evident  that 


r 


-  j'JUt 


U“(UJ)  =  Jn  Xn(r)  e  ^  dt  =  J  X~(t)  e’jUJt  dt 


i  r 


’0 


'0 


Using  Eqs .  (21)  and  (4)  for  U^a)  and  *n(a>)  ,  Eq.  (119)  becomes 


(119) 


-4=  f  x  (t)  e“ju)t  dt  =  y  \  f - 77 - rr  I 

•y2TT  Jq  m  £  mn  Ln<7  +  j(tu-n(3)j 


n=l 


(120) 


Setting  ou  =  0  in  this  last  equation  yields 

m 


r  x  (o  dt  =  ^-  y  ^ 

J0  m  CT-Jp  L  n 


n=  1 


and,  utilizing  the  identity  Eq.  (110) 

cl 

.  m 


f  x  (t)  dt  =  7^2  Miei 
m  Vtn 


m  =  1,2, .. . 


0 


(121) 


(122) 


The  important  identities  derived  in  this  section  are  summarized 


in  Table  2. 
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VI.  REPRESENTATIONS  WITH  X  (t)  AND  U  (w) 
_ m _ m 

The  preceding  sections  have  focused  on  generating  che  complete 

orthonormal  exponential  sums  X  (t)  and  Y  (t)  and  the  rational  functions 

m  m 

U  (id)  and  V  (uj)  .  Because  the  corresponding  elements  of  these  sets 
m  m 

are  related  as  Fourier  transform  pairs,  it  is  a  simple  matter  to 

determine  simultaneously  the  least-mean- squared-error  representations 

2 

of  prescribed  functions  g(t)  cL  (0,®)  by  X  (t)  or  Y  (t),  and  of 

m  m 

2 

functions  h(u))  C  L  (-»,«)  by  U  (id)  or  V  (id).  The  solution  of  the 

tn  m 

optimum  expansion  coefficients  follows  the  usual  treatment  found  in 

the  literature  on  generalized  Fourier  analysis.  Certain  noteworthy 

simplifications  evolve,  however,  because  of  the  special  nature  of  the 

orthonormal  functions  developed  in  this  Memorandum. 

2 

When  a  specified  function  g(x)  CL  (0,®)  is  to  be  approximated 

in  a  leas t-integra ted- squa red-error  sense  by  orthonormal  elements 

2 

{Q^OO]  complete  in  L  (0,®),  it  is  necessary  to  determine  the  coef¬ 
ficients  a  in  the  Mth  partial  sun; 
m 

M 

g(x)  «  £  am  0m(x)  =  g(x)  x  >  0  (123) 

m=  1 

2 

so  that  the  L  -norm 

00 

||g(x)  -  g(x) JJ 2  =  J  |g(x)  -  g(x)|~  dx  (124) 

is  minimized.  The  necessary  condition  for  achieving  an  extremum  of 
Eq.  (124)  is  that(12) 
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5(x)  -  I(x) 


!  7  y 


(125) 


or,  substituting  Eq.  (123)  fcr  g(x)  in  Eq.  (123),  that 


-  J  jg(x)  g  (x)  -  g(x)  £  aa  9a(x)  -  g*(x)  £  a_  M.x) 

I*  0  : 


M  M 


+ 1  l  %  a- 

ro=  1  n=  1 


*  * 

a  a  @  (x)  9  (x)  dx  =  0 
m  n  m  n 


r  =  1.2 . A 


Upon  differentiating  Eq.  (126)  and  utilizing  the  or t hone ma 1 i ty  of 

the  0  (x)  ,  Eq.  (126)  reduces  to 
m 


-  j  g(x)  e  (x)dx  -  j  g  (x)  3  (x.) dx  +  a_  +  a  =  0 

J0  r  -o  r  " 


(127) 


Thi  .  last  equation  in  is  obviously  satisfied  by 


r”  * 

*  •=  g(x)  a  (x)  dx 

r  r 


r  =  i}2.  . .  mM 


(126) 


The  fact  that  this  value  of  a  leadc  to  the  desired  ainicoizat ion  of 

r 

2  ( 13) 

the  norm  jjg(x)  -  g(x)jj  can  be  seen  from  the  sufficiency  ''-•ndition  1  ' 


=  12  M 
1  )*■)•••  )l* 


(129) 


For  the  particular  set  of  functions  R^(t)  given  by  Eq.  (18),  it 
is  clear  from  Eq.  (128)  that  a  least-squared-error  representation 
g(t)  of  a  function  g(t)  CL^(0,®)  of  the  form 


g(t)  *,  }  a  X  (t)  =  g(t) 
m  m 


0  <  t  <  o> 


(130) 


39 


retires  thi:  a  it  st:t  :t;  acocroir?  :c 


•r  *  r-  / —  *  ** 

g(')  X,(-)  i-  =  >  7>  V,- 


13  :> 


The  integrals  determining  these  a_  are  identifiable  as  the  La; 
ttansfjm^'^  of  g(t)  evaluated  a:  :'-t  c.nrli 
s_  =  n(r-;r),  r.  =  l,:,-.-,*,  i.e.. 


.ex  trecuenctes 


-n(o-j=)t  4. 

e  r(^.r  tt  = 


-st  .  . 
e  ?t.t  n; 


s=s 


=  c-(s  } 


•s=s 


(132) 


mo.  tnerery. 


.  w  — tr  r  -  c  - 


Consequently ,  knowledge  of  the  L-ap  lace  transtom,  C(s_)  ,  of  g(t)  at 
the  M  complex  frequencies  s_  =  n(>^_:s}  ,  n  =  1,2,...,M,  is  sufficient 
to  determine  the  expansion  coefficients  a  of  Eq.  (131. 
the  optimum  integrated-squared-error  ape roxinatior.  g(t) 
scribed  function  g(t)  on  t  >  C .  Alternatively,  this  algorithn  for 
representing  g(t)  in  terns  of  X  (t)  also  provides  a  technique  for 
obtaining  an  approximate  numerical  inverse  g(t)  to  the  prescribed 
Laplace  transform,  G(s)  ,  of  an  unknown  function  g(t). 

Through  Eqs.  ( 13C)  -  ( 132)  ,  the  approx  man  t  g(t)  of  g(t)  can  "re 


expressed  as 


M  n 


ict>.  [ 


\  P)„ 

y  ~- 


“  ~k 

K 


G(s  j  e 
nn  n 


-n(a-jp)t 


t  >  0 


(133) 


!  1  n=  1 
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As  a  consequence  of  the  Fourier  c ran* fora  relationship  between  corres¬ 
ponding  eleoents  X  (t)  and  U  («) .  a  least-integrated-squared-error 

2L  H 

-  2 

representation  h(a;  to  a  specified  function  h(x)  c  L  (-•,»)  can  also 

be  readily  computed.  If  h(«)  is  the  Fourier  transfers  of  g(t)  and  if 
h(s)  is  defined  as 

M 

h(«)  =  y  ba  » h(*>  in4> 

3>1 

with  the  coefficients  b  selected  to  vield 

21 


(fc  } 


h(«)  -  h(a) 


T  ?  w 


C35) 


then,  following  the  solution  for  the  Fourier  coefficients  a  ,  the  b 

21  Q 

are  dete rained  as 


b  « r  u*cs>  d* 

21  d  2 


a  «  1,2,. ...X 


(i3b) 


But  by  parseval's  theoren  and  the  previously  established  relation 
3(X  (t)j  »  -fir  C  (e) ,  it  ftilows  that 

£L  22 


^  J  h(i)  U  (*;  Ac 


1 

7^ 


r“  * 

g(t)  XJ  t)  dt 

•  A  ‘ — 


(137) 


sc  that  substitution  of  Eqs.  (131)  and  (135)  in  Eq.  1137)  yields  the 
numerical  equivalence 


b  *  J2n  a  m  -  1,2,. . .  ,M  (138) 

m  a 

This  result  and  another  application  of  Parseval's  theorem  car. 
provide  a  link  between  the  approximation  errors  incurred  in  the 
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dootains  (XJt)},  g(t)  CL2(0,»),  and  tU=(«)j,  h(.)  CL2(-.,.). 
c  denotes  the  error  in  the  representation  g(t)  —  g(t);  i.e. 


c  -  ain  g(t)  -  g(t) 

{•J 


2  r*  -  .2 

=  =in  j  g(t)  -  g(t)  ,  dt 

(•J  0 


( 139) 


^hen  fran  Eq.  (130)  for  g(t),  c  can  he  written  as 


r-[ 

«  *  j  j  g(t>  g 


(O  -  s  (O  Y  a=  x-(0  -  ?(=>  )  x*( t) 


M  M 


I  I  \  %  X-(°  dt 

n-l 


0*0) 


or,  taking  account  of  Eq.  (131)  for  the  oPtiau=  and  Eq.  (26)  for 
the  crthonotsality  of  the  X  (t) 

•  •  C 5  *{t>  ■**-  !•=•:-  i  «:•,*  K*: 


(i*D 


so  that 


«  =  .  g<t) 


2  r 


7  a  2 

L-  ' 

3=  1 


In  a  similar  fashion,  the  error,  £  in  the  representation  h(«) 
2 

of  h(j>)  CL  (-«,*)  can  be  derived  fros  Eqs.  (134)-(136)  and  (22)  as 


( 1*2) 


£=  sin  :  h(a’)  -  h(u>) 


h(tt)  1  -  y  i  ^  t2 


(1*3) 


Prom  the  Fourier  transform  relation  of  g(t)  to  h(o>)  ,  Parseval's 


theorem  guarantees  that 


f 


-42- 


2  1  2 
g(t)  =  —  h(«0 

Consequently,  c  and  £  arc  iclat-rJ  as 

M 

£  =  2tt  g(t)  2  -  ^  ,  aj2  =  2r  t 

3=1 

Since  j  a  !*"  is  nonnegative,  Eq.  (145)  indicates  that  both  e  and  £ 

are  aonotonically  r.onir.creat  ir.g  as  M  is  increased.  Therefore,  for 

an  arbitrarily  prescribed  approximation  error  e  or  £  M  car.  be 

iteratively  ascertained.  Moreover,  according  to  Eqs.  (131)  and  (136) 

for  a  and  b_,  M  can  be  determined  without  recomputing  any  of  the 

previously  calculated  coefficients  a_  or  b  ,  m  =  I,2,...,M-i. 

The  key  results  of  this  section  are  recapitulated  in  Table  3. 

In  the  table,  the  inner  product  notation  (x('),v(^))  denotes 
«r2  * 

x(£>  y  (;)  d",  where  the  appropriate  values  of  r  and  r  are 
»  12 

~ri 

obvious  from  the  context. 


(144) 


(145) 


REPRESENTATION  RUINATIONS 
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VII.  COMPUTATIONAL  ASPECTS  OF  ORTHONORMAL  EXPANSIONS  IN  U  (id)  and  X  (t) 

_ m _ m _ 


The  computational  labor  in  evaluating  the  coefficients  X  from 

mn 

Eq.  {45)  can  be  substantially  reduced  by  employing  a  recursive  pro¬ 
cedure.  In  order  to  develop  such  a  recursion  scheme,  the  coefficients 

A  .  must  be  expressed  in  terms  of  previously  determined  X  .* 
fefl,n  m,n 

This  can  be  accomplished  by  using  Eq.  (45). 

With  m  replaced  by  rnt-l  in  Eq.  (45),  the  coefficient  ^  becomes 

m 


I  Hgg 


mfl , 

n 


m  «=  1,2, .. . 

1  4  n  S  m  +  1 


r*l 


X  «< 

■4-1, n 


(146) 


m  *  1,2,.., 
n  >  m  +  1 


Since  X  is  alao  determined  by  Eq.  (45) ,  the  ratio  of  X  .  to  X 

mn  bh-  i ,  n  mn 

can  be  formed  as 


m  =  1 
n  =  1 


m  =  2,3,.. 
1  S  n  S  m 


(147) 


'in  order  to  avoid  confusion,  a  comma  is  used  to  separate  the 
subscripts  m  and  n  in  the  coefficients  X  .  Accordingly,  X  is 

th  mn  th  mn 

understood  to  signify  X  ,  tne  n  coefficient  of  the  m  basis 
function.  m,n 


-45- 
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and  by  employing  the  expressions 


/  oh-1  n( 


m  (rofl-n)  (a-jg)  m,n 


m  =  1,2,.. 
1  f  n  ^  m 


nH-1  ,n 


/igtilg  nstifi" 


-  (nH-1)  £ 


M-l  ,k 
k 


m  *  1,2,.. 

n  =  m  +  1 


(153) 


m  =  1,2,... 
n  >  m  +  1 


The  check-sum  relation  Eq.  (67)  still  applies  and  can  be  used  to 

detect  errors  in  the  evaluation  of  each  new  row  of  basis  coefficients 

X  .  (n  =  1,2,...,  oh-1;  m=  1,2,...)  generated  from  the  previously 
uh~  l ,  n 

computed  X 

m,n 

It  is  also  possible  to  provide  a  recursion  formula  for  the 
efficient  computation  of  the  orthonormal  elements  (X^t)}  and  {'J  (u>)}. 
From  Eqs.  (3)  and  (18),  X  (t)  can  be  expressed  as  the  sura  of 

Q 

exponentials 


X  (t) 
m 


X  e 
mn 


-n(o-jg)t 


m  — -  1,2,... 


(154) 


The  fact  that  this  sum  involves  integrally  related  decay  factors  can 

be  exploited  to  evaluate  X  (t)  recursively  for  arbitrary  values  of 

t  >  0.  This  is  accomplished  by  first  defining  the  associated  quantities 

X  as 
m,  r 


X  .  -  X 
m,  1  mm 


r  =  1 


(155) 


0 


(156) 
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and  with  X  ,  n  4  0,  given  by  the  iterative  relation  Eq.  (153).  With 
rnn 

these  definitions,  it  is  clear  that  m  iterations  of  the  expression 

x  =  e'(a'j0)t:  x  +  x  r  =  !>?’* ‘ ‘ ,ra 

m,r+l  m,r  m,m-r  m  =  1,2,... 

1 -  i-h 

and  final  multiplication  by  -jin  yields  the  m"  basis  function 


X  (t)  X  (t) 

m  m,mfl 

Thus,  for  arbitrary  t,  X  (t)  can  be  calculated  with  only  one  evalua- 

m 

tion  of  e  with  tm-1  mulLip j ieations ,  and  with  m-i  additions. 

This  is  an  important  economy  in  either  manual  or  machine  computation 
time  as  it  results  in  only  one  exponential  table  look-up  or  one 
exponential  subroutine  entry. 

Bv  following  the  calculation  of  X  (t)  ,  U  (u>)  can  also  be  itera- 

m  m 

tively  obtained  for  any  value  of  iu.  From  the  rational  function  form, 
(41),  of  U  iu)) ,  it  follows  that 


V.(a,) 


(nvf  :.)q 
TT 


m 

Jf  [no- j (w-n&) ] 

n=  1 _ 

n>f  1 

jj  [mx-j(u)-ng)] 
n=l 


=  1.2, 


Consequently, 


Sti  _ fro-  j  (ttJ-mB)  1  ,  , 

m  [(mfl)o  +  j  [  to-  (  ith-  1 )  3  3  D  m 


m  =  1,2, 


with 


UjCmO 


-/* - 1 - 

■V  rr  a  +  j  (uj-3) 


*> 


(157) 


(158) 


(159) 


(160) 


(161) 


Once  Uj(w)  is  calculated  for  an  arbitrary  value  of  w,  then  , 

follow  by  successive  multiplications  by  /2(c- juM-jfl)  /  (2o+-  Jtu-2  j  0)  , 

^3/2  (2o-jf  *2jg) /(3w-juj-3jp)  ,  and  so  forth. 

It  has  already  been  noted  that  computation  of  the  expansion 

coefficients  for  a  least- integrated-squared-error  representation  of 

a  prescribed  g(t)  requires  knowledge  of  the  Laplace  transform  G(s  ) 

‘  n 

of  g(t)  at  the  M  complex  frequencies  s^  =  n(o+j0) ,  n  =  1,2,..., M. 

Thus,  evaluation  of  the  coefficients  a  in  Eq.  (131)  entails  integrals 

m 

of  the  type 

CC 

L  {g(t)}  s  G(s  )  =  f  g(t)  e"n(t>fje)t  dt  n  ~  1,2,. ...M  (162) 

n  n  Jo 

When  the  Laplace  transform  of  g(t)  is  not  availtble,  it  can  be  approxi¬ 
mately  determined  by  a  variety  of  quadrature  schemes. 

If  the  linear  functionals  L^  are  approximated  by  the  quadrature 
formula 

M 

Ln  £  wk  g(tk )  n  =  1,2, ...  ,M  (163) 

k-1 


then  an  appropriate  criterion  for  obtaining  the  weights  wk>  k  =  1,2,...,M 


for  assigned  sample  points 

-(Ofjp)t  -2(<j+jg)t 
e  )  “  >  •  •  •  : 


t, ,t„ . t  is  that  L  be  exact  for  g(t)  =  1, 

1  /  M  n 

e  ^  This  criterion  translates 


into  the  system  of  equations 
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M 


I 


dt 


n(ofj0) 


k=l 


M 


I* 


k=l 


-(CH-jg)tk  f®  -(nfl)(0fje)t  _  _ 1 _ 

k  “  J0  (nfl)(ofjB) 


y  (164) 


x 


-(M-l)(w-jp)tk 


e-(m-M-l)  (w-jg)  tdfc 


l _ 

(ivfM-l)  (W-Jg) 


It  is  convenient  to  employ  the  following  matrix  notations  in  the 
above  equations; 


(165) 


-(CH-j  g)  t1 
e 


1 

-(W\)P)t2 

e 


V  « 


1 


-(w-jg)tM 


(166) 


-(M-D(w-jg)t 

e 


-(M-l)(ofjg)t2 


"  (M- 1)  (of  j  g)  t 


e 


e 
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n(«-j$) 


(nflKofjp) 


(167) 


- i - 

(w-m-  i)  (of  j  g) 


With  these  matr*-*-  definitions,  the  system  of  equations  in  Eq.  (164) 
can  be  written  as 


V  w  =  m 


(168) 


and  the  symbolic  solution  for  the  weights  becomes 


w  =  V  1  m 


(169) 


V  is  identifiable  as  the  Vandermonde  matrix  and  V  as  its 

inverse.*  Since  o  must  be  nonzero  and  positive,  and  since  the  sample 

points  t^,  i  »=  1,2,...,M  are  distinct,  V  is  nonsingular**  and  has  an 

inverse,  V  Once  V  *  is  evaluated,  the  solution  for  the  quadrature 

weights  can  be  obtained  from  Eq.  (169)  and  the  expansion  coefficients 

a  can  be  finally  computed, 
m 

Numerous  expressions  are  available  for  determining  the  elements 
of  the  inverse  Vandermonde  matrix.  ^  With  a  judicious  choicv' 


See  Ref.  10,  p.  92. 

» » 

See  Ref.  10,  p.  93. 
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of  the  sampling  points,  t^,  the  inverse  V  ^  can  be  computed  readily 

either  from  greatly  simplified  formulas*  or  written  directly  from 

( 19)f  * 

tables  of  Universal  matrices.  '  An  application  of  Universal 

matr'ces  to  a  quadrature  scheme  similar  to  Eq.  (163)  can  be  found  in 

Ref.  1* **  Because  of  the  direct  analogy  to  the  present  situation  (with 

a  replaced  by  (<W-j0)),  the  derivations  will  not  be  repeated  here.  By 

using  the  Universal  matrices  and  the  special  sampling  points  t.,  it 

is  possible  to  derive  the  weights  w  from  a  simple  computation  of  the 

moment  vector  m  and  a  premultiplication  by  the  known  matrix  V  ^ .  The 

expansion  coefficients  a  then  fol.ow  from  Eqs.  (163)  and  (131). 

m 

Since  a  and  b  are  related  through  Eq.  (138),  a  similar  development 
mm 

can  be  made  for  efficiently  determining  the  coefficients  b  in  the 

m 

* 

approximation  h(tu)  to  a  prescribed  function  h(u))  . 


See  Ref.  16,  pp.  96-98. 
tf 

Universal  matrices  are  the  inverses  of  the  Vandermonde  matrices 
with  sample  points  t.  =  i-n/2,  i=C,l,...,n  and  with  n  equal  to  an 

integer. 

t » t 

See  Ref.  1,  pp.  55-61. 


VIII.  SELECTION  OF  THE  COMPLEX  DECAY  CONSTANT  (a- 18) 

The  orthonormal  functions  X  (t)  and  U  (<a)  ,  the  generating  coef- 

tQ  m 

flcients  X  and  y  ,  and  the  expansion  coefficients  a  and  b  dis- 
mn  'an  mm 

cussed  in  Section  VII  are  dependent  on  a  and  0.  Except  for  the 

constraints  that  a  and  0  be  real  and  that  a  be  greater  than  zero, 

the  parameters  a  and  0  have  not  yet  been  specified. 

In  attempting  to  formulate  a  criterion  for  selecting  a  and  p, 

several  difficulties  are  immediately  encountered.  First,  in  data 

which  arise  from  a  sum  of  exponentially  damped  sinusoids  of  unknown 

decay  and  frequency  constants,  it  may  not  be  possible  to  obtain  a 

unique  complex  decay  value  (or  pole  location)  such  that  -n(o-jp), 

n  «  1,2,...,  matches  all  the  decay  factors  inherent  in  some  prescribed 

data  g(t)  ,  or  that  matches  all  the  poles  comprised  by  a  given  h(u>) . 

Second,  there  always  exist  functions  g(t)  or  h(u))  for  which  the 

optimum  choice  of  a  and  p  in  an  M-term  representation  does  not  remain 

best  as  the  approximation  complexity  is  increased. 

Another  problem  in  solving  for  the  complex  decay  constant  (a-JP) 

occurs  when  tue  given  data  is  impaired  by  noise  or  when  the  subsequent 

manipulation  of  the  data  is  accompanied  by  round-off  errors.  The 

intractability  of  this  classical  problem  is  widely  recognized  and  has 

been  amply  illustrated  even  in  ;ases  where  four  or  fewer  exponentials 

(20-30)t 

underlie  the  numerical  data.  Consequently,  the  present  objec¬ 

tive  in  solving  for  a  and  p  will  not  be  to  recover  the  original 
parameters  imbedded  in  given  data.  Instead,  an  approximation  to  the 
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data  will  be  sought  which  is  best  in  an  integral-square  sense  with 

respect  to  the  Fourier  coefficients  a  or  b  .  Another  more  tractable 

mm 

criterion  may  be  necessary  to  optimize  the  approximation  with  respect 
to  o  and  0. 

Since  a  and  0  are  subject  to  only  the  cwo  aforementioned  con¬ 
straints,  their  choice  is  essentially  arbitrary  and  can  be  based  on 
any  criterion  of  optimality.  One  obvious  criterion  is  the  minimisation 
over  a  and  0  of  the  integral- square  approximation  error  t  or  ^  given 
by  Eqs.  (142)  and  (14_>).  Although  this  objective  would  be  consistent 

with  the  conditions  lead  ,.g  ft  the  expansion  coefficients  a  and  b  , 

mm 

it  unfortunately  results  in  a  nonlinear  programming  problem  requiring 
iterative  search  procedures  for  its  solution.* 

A  more  tractable  criterion  than  least-squares  is  one  which 
requires  matching  the  asymptotic  approach  to  zero  of  the  approximant 
and  prescribed  function  for  large  t.  More  precisely,  if  g(t)  and  itB 
derivative  g'(t)  exist  for  some  t  »  1  and  are  nonzero,  both  ct  and  0 

A 

can  be  determined  by  requiring  g(t)  and  its  approximant  g(t)  to  have 
the  same  decay  envelope  for  large  t.  This  condition  can  be  derived 
from  Eqs.  (130),  (18),  and  (3),  as  follows. 

Since 

M  Mm 

sco  „  g(t) .  y  .  x  <t> .  y  ,  y 

m  m  LmL  mn 

1  m=  1  n*  1 


t 

See  Ref.  1,  p.  40. 


t  >  0  (170) 
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the  first  derivative  of  g(t)  can  be  written  as 

M  m 

8'(0  a  «'(t)  -  -  -v/*T  (o-Jg)  Y  a  \  n\  e‘n^CT'JP)t  t  >0  (171) 

Zj  m  4j  mn 

n>“l  n*l 

For  sufficiently  large  t  and  with  a  >  0,  it  is  clear  that  only  the 
terms  with  the  smallest  decay  factor,  -a,  predominate  in  Eqs.  (170) 
and  (171).  Consequently,  Eqs.  (170)  and  (171)  can  be  simplified  for 
t  »  1  to 

M 

g(t)  m  g(t)  *V2tt  e“(0"jp)t  7  a  t»l  (172) 

[_ j  m  ml 

m*l 

and 

M 

g'(t)  m  g'(t)  «  (o-jg)  e'(CT'jP)t  y  a  \  t»l  (173) 

/ .  m  mi 

m=l 

i.f  the  ratio  of  these  last  two  equations  is  formed,  then  for  a  >  0 


g(t) 

"  !(t) 

mt  o-jg 

g(t) 

t»i* 0 

(174) 

t»i 

t»l 

Hence , 


a  mt  -Re 


ft'(t) 

g(t) 

t»l 

g(t)  4  0,  g'(t)  f  0 


(175a) 


and 


P  w  Im 


J 

g(t) 

t»l 

(175b) 


When  the  prescribed  function  g(t)  is  real,  the  solution  for  g 
given  by  Eq.  (175b)  indicates  that  g  can  be  set  equal  to  zero.  However, 
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this  may  be  an  unsatisfactory  choice  for  3  if  g(t)  exhibits  an  oscil¬ 
latory  nature  for  some  values  of  t.  In  such  cases,  3  and  o  can  be 
resolved  by  representing  g(t)  not  as  in  Eq.  (170),  but  as 
M  M 

g(t)  *  g(t)  =  T  c  Y  (t)  =  V  c  [X  (t)  +  X*(-t)] 
m  m  //  m  m  m 

m= 1  m=  1 


where  the  orthonormal  functions  Y  (t)  are  defined  by  Eqs.  (55),  (61) 

m 

A 

and  (62).  In  view  of  Eq.  (62),  g(t)  becomes 
M  m 


g(t)«^  l  cm  £  |xj 


-na 


m=-l  n=l 


l^^cosf^t  i-  -j~y 


+  -[f|  j  sin[n3|t|  +  arg(Xmn)] 


arg(>J 


t  <  CD 


For  t  >'  1,  the  terms  for  which  n  =  1  predominate,  as  before,  so  that 

N 

I 

cos  [3t  +  arg(\ml)] 


-at 


8(t)„/„e~  £  g*mlM 


with 


m=l 


+  j  3in(jt  q-  arg(\ml)] 


M 


t»l 


g'(t)  «  -a  g(t)  +  /n  e‘at  Y  cmUml  I 


m=  1 


+  j0  cosfpt  +  arg(\ml)]l- 


-3  sin[3t  +  argUml)] 


t»l 


(176) 


(177) 


(178) 


(179) 


and 


-56- 


M 


g'(0%-0g'(t)+,/iTe'Otffp  £  cml\nlh  8lnC3t  +  arg^ml)3 


M 


-  i  co.[et+  .rgaml)]U  A  82  £  cjxj 


m»l 


cos[pt  +  arg(Xml)] 


+  J  sin[0t  +  arg(Xml)] 


t»l 


(180) 


or 


g"(t)  *  -2og'(t)  -  (o2+32)  g(t) 


t>>l 


(181) 


Similarly, 


g*(t)  «  -  2c  g"vt)  -  (a2fg2)  g'(t)  t»l 


(182) 


The  preceding  two  equations  relate  the  unknowns  o  and  8  to  the  given 

function  g(t)  and  its  derivatives  at  some  large  value  of  t.  Equation 

2 

(182)  can  be  solved  for  0  to  yield 


.  slSi 1  +  2qg"(t)  +  cr2z'(t) 

g'OO 


t»l 


8  /(t)  |t»l  ^  °  (183) 


2 

When  this  expression  for  0  is  substituted  into  Eq.  (181),  the  follow¬ 
ing  relation  for  a  is  obtained 


g"(t)g(t)-g'(t)g*(t) 
2[g  ,2(t)-g(t)g"(t)  ] 


g'(t)  ^±->/8(t)g'(t),  t»l  (184) 


Thus,  0  is  also  explicit  iy 


related  to  g(t)  and  its  derivatives  for 


some  large  value  of  t  as 
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0  *  V 


S^l  .  SlL*±  . 

g,(t)  8  (t)  g'2(t)  -  g(t)g"(t) 

4[g/4(t)+g2(t)g"2(t)-2g'2(t)g(t)g/,(t)] 


(165 


t»l 


g'(t)  /  0,  g'(t)  ^  ±R(c)g'(t), 


t»l 


Several  qualifications  are  noteworthy  regarding  these  equations 

for  ct  and  8-  It  is  clear  that  Eqs.  (175),  or  Eqs.  (184)  and  (185), 

provide  suitable  values  of  (a-jg)  when  the  prescribed  g(t)  i3  given 

as  a  long-time  record.  When  g(t)  is  a  transient  or  pulse-type  function, 

or  when  the  complete  time  history  of  g(t)  is  unknown,  tie  asymptotic 

properties  of  g(t)  are  not  available  for  estimating  ct  and  8  by  the 

aforementioned  equations.  For  these  pulsatile  functions,  it  is 

necessary  to  resort  to  more  elaborate  techniques  for  obtaining  (ct-J8) . 

One  procedure  for  obtaining  a  and  g  relies  on  the  discrete  version 

of  Prony's  exponential  approximation  method.1  In  essence,  an  applica- 

(31) 

tion  of  the  J  ony  algorithm  enables  a  selection  of  (ct-J8)  based  on 

many  samples  of  g(t),  rather  than  on  a  match  of  only  the  asymptotic 

values  of  g(t)  and  its  derivatives  for  large  t.  Other  versions  of 
(32-33) 

Prony's  scheme  which  are  available  can  also  be  used  to  select 

A 

ct  and  8.  Depending  on  the  procedure  chosen,  the  approximant  g(t) 
satisfies  in  a  least-squares  sense  a  finite  difference  or  differential 
equation  involving  the  sampled  or  continuous  ordinates  of  the  specified 
g(t).  Since  all  these  approaches  are  conceptually  similar,  and  since 

See  Ref.  1.  pp.  42-51. 
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the  discrete  version  of  Prony's  method  has  been  detailed  in  Ref.  1 

for  the  case  in  which  g  ■  0,  the  derivations  will  not  be  repeated 

here  with  g  ji  0.  It  will  suffice  to  state  that  the  derivations  in 

Ref.  1  apply  to  the  condition  g  ft  0  considered  in  this  Memorandum 

by  merely  substituting  (o-jg)  for  a  in  Eqs.  (127) - ( 154)  of  Ref.  1. 

Though  no  mention  has  yet  been  made  of  the  associated  problem 

of  selecting  a  and  g  for  approximations  of  functions  h(<u)  c  L  (-»,») 

by  sums  of  the  orthonormal  elements  U  (u>)  ,  it  turns  out  that  all  the 

m 

techniques  discussed  so  far  are  applicable,  albeit  indirectly  and 

A 

with  additional  computational  labor.  Since  the  approximants  h(u>) 

A 

and  g(t)  are  Fourier  transform  pairs,  the  Fourier  transform  can  be 
taken  of  a  graphically  or  analytically  specified  h(u>)  to  produce  a 

g(t)  or  samples  of  g(t)  at  S  points  t  «  0,1 . S- 1  .*  In  turn,  these 

ordinates  of  g(t)  can  be  utilized  in  all  of  the  Prony  schemes  cited 
earlier,  and  a  choice  of  (a-jg)  can  again  be  made  based  on  g(t)'s 
satisfying  a  difference  or  differential  equation. 

See  Ref.  30,  pp.  67-75. 
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IX.  NUMERICAL  EXAMPLES  AND  APPLICATIONS  TO  FILTER  DESIGN 


In  the  introduction,  several  important  practical  applications 

are  enumerated  for  the  orthonormal  sets  X  (t),  Y  (t),  U  (at),  and 

m  m  m 

V  (oj)  .  In  order  to  clarify  the  use  of  the  algorithms  developed  in 
m 

this  Memorandum,  two  filter  design  problems  will  be  illustrated  in 
this  section. 

The  first  example  pertains  to  an  approximation  problem  in  network 
(34) 

synthesis.  It  requires  finding  a  physically  realizable  transfer 

function  for  a  finite,  lumped-element,  passive,  linear  network  such 

that  the  network's  impulse  response  is  a  replica  of  the  waveform 

depicted  in  Fig.  1.  The  approximating  response  must  be  close  enough 

to  the  prescribed  response  that  the  mean- square  error  between  them 
_2 

is  less  than  £  X  10  over  the  time  interval  (0,10).  In  quantitative 
^rms,  the  prescribed  transient  resronse  is  given  as 


0 


2t 


g(t) 


-2(t-l) 


-2.42377t 

e 


t  <  0 

0  S  t  s  0.5 
0.5  St*  0.95 
t  >  0.95 


(186) 


A 

and  it  is  necessary  to  find  and  N  for  the  approximart  g(t)  such 
that 


k*  i  k«  1  1 


(187) 


and 


If10  .2 

«N  -  To  J0  t  dt  *  6  -  « 


x  10 


-2 


(188) 


Graph  of  preicrlbed  g  ( t ) 
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In  order  to  meet  this  design  specification,  the  coefficients  a^, 

k  «  1,2,..., N  oust  be  computed  for  some  choice  of  o,  3>  and  N.  In  so 

doing  it  is  convenient  to  recall  that  the  expansion  coefficients  a^ 

are  independent  of  N.  Thus,  N  can  be  selected  arbitrarily,  cu  can  be 

evaluated  and  compared  with  6,  and  N  can  be  iteratively  adjusted  to 

be  smaller  or  larger  depending  on  whether  c  <  6  or  t  >  6.  In  either 

N  n 

cate,  the  a^  are  determined  from  Eq.  (131)  aa 

k 

\  «  J  g(t)  X*( t)  dt  «  y  f  g(t)  i*a(afJ3)t  dt  (189) 

^  0*1 

k  =  1,2, ... ,N 


Thus,  evaluation  of  all  he  a^  entails  integrating  the  N  quantities 


g  -  r  dt 

■a  JQ 


m  *  1 ,2 , . . .  ,N 


(190) 


i.t 


s  =  2  r 

<E  J  , 


.0.5  ,  .0.95 

t  dt+ ,f  (1.t) 

J0.5 


e-ai(Of  j3)  t.  dt 


(191) 


j: 


e -(nCH-mj 3^-2 .42377)  t 


0.95 


2  2 
<“  (W-J3) 


_4e-0.5«(«.JS)  +  f2  .  .1(ff+jW>]  ,-0-*SmW» 


>  ♦  21 


(192) 


m  =  1,2 , . . . ,N 


In  order  to  determine  the  moments  G  (the  sampled  Laplace  trano- 

tn 

form  of  g(t))  by  this  last  relation,  a  selection  of  a  and  3  is 
necessary.  Since  g(t)  is  a  pulse-like  function  with  an  exponential 
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tail  and  no  oscillatory  character,  Eq.  ^  17 S)  is  appropriate  for  deter¬ 
mining  a  and  g.  Accordingly 


a  m  -Re 


• 

r 

" 

Alt!  j 

=  -  Re 

-2.42377e‘2,42377t 

g(t)  |t»l 

-2. 42377 t 

t»l 

_  .. 

.  C 

• 

(193) 


or 


CT  «  2.42377  >  0 


(194) 


and 


8  M  Im 


£. 'AH 

■ 

g(t) 

t»l 

r,  0 


l 

With  these  choices  for  a  and  g ,  the  first  nine  coefficients  a^  are 
found  from  Eqf.  (192),  (189),  and  the  tabulated  \  of  Appendix  A 

as 


ax  =  0.3738  +  j  0. 

a2  =  0.3957  +  j  0. 

a3  =  0.0694  +  j  0. 

a  =  -0 . 1164  +  j  0 . 
4 


a5  =  -0.1227  +  j  0. 
a  =  -0  0412  +  j  0 . 

D 

a?  =  0.0284  +  j  0. 

ag  =  0.0441  +  j  0. 


a9  =  C .0227  +  j  C. 

These  values  for  can  be  inserted  in  Eq,  (187)  along  with  the 

A 

tabulated  \  to  find  g(t)  for  t  >  0.  When  this  is  done,  the  rms 
mn 

error  over  the  interval  (0,10)  can  be  computed  from  Eq.  (188)  to 
give 


(195) 


(196) 


*9  =  0.2  x  10  =6 


N  =  9 


(197) 


The  values  of  listed  in  Eq.  (i96)  have  been  roune’ed  to  four 
significant  digits. 
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-2 

Since  60  =  0.3  X  10  >6,  it  is  clear  that  the  solution  for  N  =  9, 

O 

a  =  2.42377,  p=0,  satisfies  the  design  specification.  For  these 

* 

values  of  N,  a,  and  g,  graphs  of  g(t)  and  g(t)  can  be  compared  (see 
Figs.  2  and  3) . 

Finally,  since  the  Laplace  transform  of  X^(t)  is 

m 

r”  *  (t>  e~st  dt .  y  &  x  r  e<-no-M-jnWt  dt 

JQ  m  L  mn 

n=l 


m 


mn  sfn(a-jg) 


- 

s  =  or  +  j  g 
■s  a  +  na  >  0 
n  =  1,2,... 


(I98a) 


the  approximating  realizable  network  transfer  function  can  be  written 
as 


N  m 

ws>  -  l  l  *>  H's)  *  f 8(t)  e'st  dt  <198b> 

m=  1  n=l  ° 


Thus,  with  a  equal  to  a  real  number  and  g  equal  to  zero  (as  in  the 

A 

present  example)  H(s) ,  the  network  approximant,  can  be  implemented 

with  either  RC  or  RL  elements  to  give  the  impulse  response  g(t)  g(t). 

The  second  numerical  example  deals  with  a  problem  of  optimal 

filtering.  In  processing  signals  impaired  by  additive  noise,  it  is 

possible  to  accomplish  smoothing  and  prediction  by  constructing  an 

(35) 

appropriate  filter.  Usually  this  requires  approximating  the 

spectral  densities  of  the  signal  and  random  noise  process  by  rational 
functions.  This  procedure  in  turn  permits  the  analytically  optimum 
filter  to  be  approximated  in  the  form  of  a  lineal.  lumped-constant 
network.  Since  the  power  spectral  density  is  a  rational  function  of 
frequency  for  a  signal  whose  correlation  function  is  a  sum  of 


(4)6  junuuyojddD  S4|  puD  (4)6  u»a|£) 


Fig .  2  —  Graphs  of  g  ( t )  and  g  ( t )  for 


Re  [  g  ( t )  ] 


Fig  .3  —  Graphs  of  g  ( t )  and  g  ( t )  for  N 
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exponentials,  it  is  possible  to  use  advantageously  the  sets  {X^(t)j 

and  {Y^(t}}  in  the  approximation  of  empirical  correlation  data. 

A  commonly  observed  autocorrelation  function  for  a  stationary 

/  *1£\ 

random  prwt  <;  s  the  exponentially  damped  cosine''  1 


w/  \  -2  T 

Y(t)  =  e  11  cos  nr 


T  |  >  0 


(199) 


Though  the  best  approximant  of  Y(t)  of  exponenticl  form  is  clearly 
Y(t)  itself,  it  is  instructive  to  see  how  efficiently  this  exponen¬ 
tially  damped  cosinusoid  can  be  approximated  in  an  integral-square 

sense  bv  the  orthonormal  elements  X  (t) .  Since  Y(t)  is  defined  over 

m 

the  doubly  infinite  interval,  I T i  <  ®,  and  since  the  X  (t)  are  non- 

*  (Q 

zero  only  over  t  >  0,  it  is  necessary  to  represent  Y(t)  over  (-oo,a>) 


f(T>  -  f+(T)  -  l  ak  X.  (T) 


0  £  T  <  oo 


(200) 


Y(t)  as  Y_(t)  =  Y  \  Xk(T) 


<  T  S  0 


(201) 


It  is  evident  from  Eq.  (199)  that  Y(t)  is  an  even  function  of  T; 
i.e.  , 


Y(t)  =  Y(-t) 


(202) 


Hence 


Y+(t)  =  Y_(-t) 


(203) 


ak  3  bk 


k  =  1.2 . N 


(204) 
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Consequently ,  fhe  moments  G  for  Y(t)  can  be  laced  simultaneously 

m 

to  a  and  b  through  Eqs.  (189),  (190),  and  (204)  as 

K  K 


G 

m 


e-(mcH-jm8f2)T 


cos  TTT  dT 


tn(w-IB)  -f  2 
[m(cH-jg)+2]2  +  T72 


m  =  1,2, ...  ,N 


In  terms  of  these  G  ,  the  expansion  coefficients  for  Y  (t)  and 
m  + 

A 

Y  (t)  become 

k 

a,  =  Jin  Y  \*  G  =  b,  k  =  1,2,. ..  ,N 

k  v  km  m  k 

m=l 

* 

where  the  are  tabulated  in  Appendix  A  as  rational  functions  of 
a  and  3. 

A 

In  order  to  complete  the  approximation  Y(t)  em  Y(t) ,  a  selection 
of  the  parameters  a  and  g  is  necessary.  This  will  enable  the  evalua¬ 
tion  of  G  ,  a.  ,  and  b  of  Eqs.  (205)  and  (206)  and,  thereby,  of 
m  k  k 

Y+(t)  and  Y  (t)  in  Eqs.  (200)  and  (201).  Since  Y(t)  is  given  as 
a  real,  oscillatory,  exponentially  decaying  function,  Eqs.  (183)-(185) 
can  be  used  to  obtain  a  and  g,  with  g  not  necessarily  zero  (as  would 
be  the  case  if  Eq.  (175)  were  used).  Accordingly,  the  first  three 
derivatives  of  Y(t)  for  t  >  0  become 

Y'(t)  =  -e  2t  [2  cos  ttt  +  n  sin  ttt] 

-  2t  2 

Y"(t)  =  e  [(4-n  )  cos  ttt  +  4tt  sin  ttt] 

Y"'(t)  =  e  2T  [(-8f6n2)  cos  ttt  +  (tt2-12)  tt  sin  ttt] 


(205) 


(206) 


(207) 

(208) 


(209) 
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(213) 

For  f (t)  prescribed  in  Eq.  (199),  the  values  am  2  and  0  >*  n 

just  derived  have  special  intuitive  appeal.  When  these  values  are 

used,  and  when  a  mean-square-error  tolerance,  Eq.  (188),  of 
-2 

6  »  0.9  x  10  is  selected,  it  is  necessary  to  compute  twelve  expan¬ 
sion  coefficients.  Thus,  with  N  -  12,  a  »  2,  and  0  -  n,  the  of 
Eq.  (189)  become 

k 

*k  '  J*”  l  \m  G«  '  bk 

n»  1 


k  -  1,2,. ...N 


(214) 
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where  now 


G  .  f  T(T)  e-m(0+J»)T  dT  -  r  c0,  „  dT 

®  » n  •  n 


or 


3  . _ aigtLBLt  2, 

m  [m(«-jg)+2]2  +  tt2 


m  *  1 1 2 1 ■ ■ i  |N 


(215) 


(216) 


Thus,  the  expansion  coefficients,  rounded  to  four  significant  digits/ 
are 


3 

S4 


0.3221  -  J  0.1132 
-0.0967  +  j  0.0452 
0.0798  +  j  0.0281 
-0.0242  -  j  0.0648 


a5  =  -0.0316  +  j  0.0488 

a£  -  0.0500  -  j  0.0014 

o 

a7  =  -0.0249  -  j  0.0361 
ag  ~  -0.0154  +  j  0.0358 


ag  =  0.0347  -  j  0.0048 

a  =  -0.0197  -  j  0.0249 
10  (217) 

an  =  -0.0107  +  j  0.0270 

al2  =  0.0265  -  j  0.0037 


with 


12 


0.895  X 


10'2  <  6 


(218) 


and 

*u  =  0.934  X  10'2  >  6 


(219) 


The  corresponding  approximation  ¥(t)  and  the  given  autocorrelation 
function  Y(t)  are  compared  in  Figs.  4-7. 

/  i  ^  \ 

Finally,  since  the  Wiener-Khintchine  Theorem'  '  relates  the 
autocorrelation  function  and  spectral  density  as  Fourier  transform 
pairs,  the  power  density  spectrum  associated  with  the  approximation 
12 

VT)  “  I  \  VT)  *  f(T)  T  2  0  (220) 

k-1 

. 

Sixteen  significant  digits  were  computed  in  all  the  calculations 

of  a,  since  the  \  may  become  large  for  certain  values  of  o,  0,  and  m. 
k  mn 


IT) 

CM 


CM 


II 

z 


CM 

II 

z 


■'f 

O) 


Graphs  of  the  real  parts  of  'fr(r)  and  'i'(r)  for 


rid  ^(t)  for 


Fig. 6 —  Graphs  of  the  real  part s  of  'i'(r)  and  'i'(r)  for 


)  and  ^(r)  for  N 
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?_(t)  “  ^  ak  VT>  -  *(t) 


T  S  0 


(221) 


*(u>)  =  J  Y+(t)  e‘JU,T  dr  +  J  J.(T)  e'ju,T  dT 


(222) 


pa  J  't'(T)  e  •^UJT  dT  =  i(w) 


Thus,  by  changing  the  variable  in  the  second  integral  of  Eq.  (222)  and 
by  noting  Eq.  (221) ,  Eq.  (222)  becomes 


*(u>)  £  ak  {J  \(T)  e'ja)T  dT  +  J  Xk(T)  eja,T  dTj 


(223) 


Since  X  (t)  and  U  (u>)  are  related  as  Fourier  tr-r.sforms  (to  within 

K  K 

a  scale  factory  see  Eq.  (23))  ,  Eq.  (223)  ran  be  written  as 


i(u>)  «  Y  ■yin  a  [U  (u>)  +  U  (a1)] 
km  m 


(224) 


or,  in  view  of  Eqs.  (21)  and  (4),  as 


12  k 


k=l  n=l 


, _  -  ~  no  Re[l  ]  +  (u>-n0)  Im[>  ] 

*(uj>  m  2  >  y  a  < - - ~ - -j - —>«#(»»>) 

L  L  k>  (no)  +  (w-n3)  ' 


Thus,  for  any  value  of  N,  once  the  a^  are  determined  for  the  auto¬ 
correlation  function  Y(t) ,  the  approximating  spectral  density  can 
be  written  by  inspection.  In  view  of  Eq.  (143)  ,  moreover,  the 

A 

integral-square  error  between  4(u>)  and  4(iu)  in  Eq.  (225)  is  simply 


(225) 


2nv 
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Appendix  A 

THE  ORTHONORMAL  BASIS  COEFFICIENTS  \ 

mn 


The  rational  function  of  a  and  3  given  in  Eq.  (45)  can  be  used 

to  generate  the  coefficients  \  .  As  indicated  in  Section  VII,  how- 

mn 

ever,  there  is  considerable  computational  advantage  in  recursively 

evaluating  \  by  Eq.  (153).  A  program  for  doing  this  is  presented 
mn 

in  Appendix  B. 

In  order  to  tabulate  the  algebraic  \  obtained  by  computer,  the 

mn 

following  definitions  are  convenient.  From  Eq.  (45),  \  is  expressed 

mn 

as  the  proper  rational  function. 


\  -  -S 

mn 


m- 1 


jj  [n(a-jp)  +  r(of j  3) ] 


(o-j0)m  1  jj  (r-n) 
r=  1 


m 


If  T]  and  A  are  defined  as 
1  mn 


n  =  m  =  1 


m  =  2,3,... 

1  £  n  £  m  (226) 


m  =  1,2,... 

n  >  m 


T1  *  j|3 


and 


(227) 


m-  1 


m-  i 


Jj  [n(a-J3)  +  t(w-)3)j  JJ  [n(a-Tj)  +  r(afT.)] 


mn 


Ul  / 

ft  <r-") 


[j  ( r-n) 


(a-Tj)m' 1  a  (228) 
mn 


r-1 


r»  1 
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where  cr  Is  given  by  Eq.  (68)  ,  then  \  can  be  reexpressed  In  terms 
mn  mn 

of  the  polynomial  A  in  o  and  as 
mn  1 


n  =  m  =  1 


m  =  2,3,... 
1  S  n  S  m 


(229) 


m  =  1,2, , 
n  >  m 


The  quantities  A  defined  by  Eq.  (228)  are  shown  in  Table  4 
mn  ' 

for  the  range  m  *  1,2,..,, 10  and  n  -  l,2,...,m.  For  each  value  of  m, 

the  check-sum  relation  Eq.  (67)  holds,  thereby  verifying  the  exactness 

of  the  rational  \  derived  from  the  polynomials  A  . 

mn  mn 

The  A  of  Table  4  have  been  generated  by  encoding  the  recursive 
mn 

relation,  Eq.  (153),  in  ALTRAN,  a  symbolic  manipulation  language,  and 
executing  the  program  o^  a  digital  computer. The  ALTRAN  compiler, 
which  i s  required  for  execution  of  the  program  listed  in  Appendix  B, 
produces  MAP  output  consisting  of  transfers  to  ALFAK,  a  group  of  sub¬ 
routines  for  computing  and  simplifying  polynomials  and  certain  rational 
functions . 

Since  ALTRAN  is  not  widely  available  and  since  it  requires  large 

computer  storage  for  evaluating  A  for  m  as  high  as  ten,  another 

mn 

program  is  given  in  Appendix  C  for  computing  the  A  .  The  Fortran  IV 

mn 

routine  listed  in  Appendix  C  computes  the  integer  coefficients  of  the 
variables  9  and  0  in  the  polynomials  A  .  The  computational  basis 
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for  the  program  is  discussed  in  Appendix  D  where  a  series  representa¬ 
tion  for  A  is  derived  to  replace  the  product  form  given  by  Eq.  (228). 
mn 

Finally,  in  Appendix  E  a  Fortran  IV  program  is  presented  which  enables 

2 

arbitrary  functions  g(t)  CL'(0,“)  to  be  expanded  in  terms  of  the 
orthonormal  elements  X  (t) . 
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Table  4 

THE  ASSOCIATED  ORTHONORMAL  BASIS  COEFFICIENTS  A 

mn 


A21  =  2a 

1^22  =  “3cr  +  T) 

2 

A31  =  3a  +  aT| 

2 

A^  =  -12  a  +  4ai] 

2  2 
A^  =  10a  -  7oP,  +  7] 

A^  =  (12ct3  +  10a2  71+  2o7]2)/3 

3  2  2 
A42  =  -30a  +  4ct  7]  +  2a7) 

A43  =  60a3  -  42a2 7]  +  6aT)2 

A4A  =  (-105a3  T  113a271  -  35a712  +  3713)/3 

A^  =  (30a4  +  43a3  7]  +  20o2  7}  +  3a7f)/6 

^  3  2  2  3 

A52  =  (-180a  -  36a  7]+  20a  71  +  4a71  )  /3 

A53  =  210a4  -  117a371  +  W 

A^  =  (-840a4  +  904a  71  -  280a  71  +  24a71  )/3 

A^5  =  (756 a4  -  110 la3 7)  +  536a2712  -  lOlaT}3  +  6T)4)  /6 

A61  =  (90a3  +  189a  *7}  +  146a  71  +  49a  Tf  +  6<771  )/15 

A62  =  (-315a5  -  198a471+  8a3 712  +  22a2 713  +  301J4)  /3 

5  4  3  2  2  3  4 

Ag3  =  560a  -  172a  71  -  78a  7)  +  8a  71  +  2aTl 
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Table  4--continued 

=  (-3780a5  +  3648a4 T)  -  808a3T|2  -  32a2Tj3  +  12aTl^) /3 
A^5  =  (3780a5  -  5505a4 T)  +  2630a3 T)2  -  505a2T]3  +  30aT)4)/3 
A66  =  (-2310a5  +  4247a4 T]  -  2842a3?2  +  852a  V  -  ll2aT)4  +  5Tl5)/5 

C  c  /O  0/  C 

A?1  =  (630a  +  1773a  T|+  1967a  T]  +  1073a  T]  +  287 a  T\  +  30aT)  )/90 
A72  =  (-2520 a6  -  2844a5 T]  -  728a4  !]2  +  208a3 T]3  +  112a2!]4  +  12a!]5) /15 
A73  =  (2520a6  +  66a57]  -  609a4T|2  -  81a3 7]3  +  21a2T]4  +  3oT)5)/2 

6  S  a  2  3  3  2  A  q 

A74  -  (-37800a  +  28920a  T1  -  784a  T]  -  1936a  ?,  +  56a  T]  +  24a!]  ) /9 

A75  =  (41580a6  -  56775a5!]  +  23975 aV  -  2875a3!]3  -  175a2!]4  +  30aT]5)/6 

A?6  =  (-27720a6  +  50964a5!]  -  34104a4!]2  +  10224a3!]3  -  1344a2!]4 
+  60a!]5) /5 

A?7  =  (154440a6  -  343 146a5  T]  +  293243a4!]2  -  122023a3 T|3  +  25703a2!]4 
- 2547a!]5  +  90T]6)  /90 

AgL  =  (2520 a7  +  8982a6!]  +  13187aV  +  10193a4!]3  +  4367a3!]4  +  OSUM]5 

+  90a!]6) /3 15 

Aq2  =  (-11340a7  -  19098a6!]  -  10386a5!]2  -  884a4!]3  +  1024a3!]4  +  334a2!]5 
+  30a!]6)/45 

7  6  S  2  a  3  ^  /  25 

Ag3  =  (12600a  +  5370a  T]  -  2913a'!]  -  1623a  !]  -  57a  T]  +  57a  T) 

+  6a!]6) /5 

Ag4  =  (-103950a7  +  51180a6!]  +  19534a5?2  -  5912a4?3  -  1298 aV 
+  108a2!]5  +  18aT|6)  /9 
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Table  4--continued 


Ag5  -  (249480 o7  -  299070c6]]  +  87075o5!|2  +  6725oV  -  3925a3  !]4  +  5oV 
+  30cT]6)/9 

Ag6  .  (-180180a7  +  317406a6!]  -  196194a5!]2  +  49404o4T)3  -  3624a3!]4 

-  282a2!]5  +  30oT]6)  /5 

Ag7  =  (1081080a'  -  2402022a6!]  +  2052701aV  -  854161a4!]3  +  179921a3!)4 
-17829a2!]5  +  630a!]6) /45 

Ag8  =  (-2027025 o'  +  5282325a6T(  -  5503581a5Tj2  +  2947489a4!]3  -  867299a3]]4 
+  1383 19aV  -  108b 3a!]6  +  315]]7)/315 

A^  =  (22680 a8  +  98478a7!]  +  181557a6]]2  +  184046a5T]3  +  110654a4]]4 
+  39398a3!]5  +  7677a2!]6  +  630CT]7) /2520 

Ag 2  =  (-113400a8  -  259020a7!]  -  218448a6 !]2  -  71156a5 T]3  +  4936a4]]4 
+  9484a3]]5  +  2304a2!]6  +  180a!]7) /315 

Ag3  =  (46200a8  +  40690a7!]  -  1731a6!]2  -  10806a5]]3  -  2914a4 ]f  +  114a3 T|5 
+  117a2  T]6  +  10a!]7) /10 

Q  A*?  c  // 

Ag^  =  (-1247400a  +  1983600  T]  +  439  128a  !]  +  7192a  1)  -  39224a  J\ 

-  3896a3]]5  +  648a2]]6  +  72oT]7)/45 

Ag5  =  (2243240a8  -  3139470G7]]  +  234765a6 T]2  +  348650a5]]3  -  30850a4]]4 

-  11710a3!]5  +  405c2  Tj6  +  90a!]7) /36 

Ag6  =  (-840840a8  +  1361108a7]]  -  703968a6 if  +  99756a5]]3  +  16024a4 T]4 

-  3732a3!]5  -  48a2!]6  +  20a!]7) /5 

-  (16216200a8  -  34949 2 50a7 T.  +  2"388493o6T]2  -  10759 7 14cV 
+  1844654a4 T]4  -  87514a3]]5  *  8379a2]]6  +  630cT]7)/90 
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Table  4--continued 


=  (-32432400*?  +  84517200o/T)  -  88057296c°7)  +  47 159824(7  Tf 

-  13876 7 f  t«j* T]4  +  2213104c3!]5  -  1 73808a2 !]6  +  5040on?)  /3 15 

*  (6806800a8  -  20355850a77|  +  25047613c6  T]2  -  16458298a5Tl3 

//  C  /•  ^  Q 

+  6267166a  7)  -  1402306a  T)  +  177733a  T)  -  11458011  +  2801)  )/l 
=  (11340a9  +  5831 10a8  T]  +  130 1697a7T)2  +  1646458a6  T|3  +  1289454o5l/ 

Ac  ^  /■  0  0  O 

+  639606a  T\  +  195977a  T)  +  33858a  1)  +  252001]  )/11340 

=  (-311850a9  -  910755a8Tj  -  1054017a77)2  -  577963a6 T]3  -  110949a5l/ 
+  34719a4 1]5  +  22933a3!]6  +  4527a2!]7  +  315al]8)/630 

=  (277200a9  +  382740a8!]  +  111684a7 1]2  -  70029a6!]3  -  49902a5T|4 

-  8058a4!]5  +  1044a3!]6  +  411a2 T,7  +  30aT]8) /35 

-  (-8108100a9  -  1829160a8Tj+  3350232a7Tj2  +  1144568a6!]3 

-  236976a5!]4  -  123384a4!]5  -  5528a3 T)6  +  2088a2!]7  +  I80a!]8)/135 

=  (4540536a9  -  3097962a8!]  -  927117a7!]2  +  582016a6!)3  +  96270a5T4 

-  28734a4!]5  -  4117a3!]6  +  288c2T,7  +  J6af)8)/18 

-  (-3153150a9  +  4473525a8!]  -  1619049a7!]2  -  153891a6!]3 

+  134907CV  -  1977a'  !]5  *  2979a3T6  +  39c  V  +  15oT)8)/5 


A10  7  *  ( 1297 29600a9  -  263377800a8!]  +  192  158694a7!]2  -  57689219a6!]3 

+  39975  18a5!]4  +  1144542a4!]5  -  154546c3T6  -  3339aV  +  630a!]8)  / 135 

A10  8  *  *  2756  75400a9  +  702 180000a8ri  -  706228416a7!]2  +  356829856a6!]3 
-  943  72752a5!]4  +  11872992a*!]5  -  370816a3’-6  -  44064c2T7 
+  2520aT]8)  /  J  1 5 
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Table  4--cont inued 

A10  9  “  (61261200°9  -  183202650o8T1  +  225428517oV 
+  56404494c5  Tj4  -  12620754c4 J\5  +  1599597c3T]6 
+  2520CT]8)  / 140 

A10  10  =  ("104756652c9  +  353598543c8 T)  -  502107309c7T1 

-  18384026  1o5T]4  +  5344374 3c4T5  -  9532127c3 T\ 

-  53955c T]8  +  1134 T]9)  / 1 134 


-  148124682c6  Tj3 
-  103l22cV 

2  +  39 16  7286  5c8  T)3 
6  +  993411o2Tj7 
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Appendix  B 

AN  ALTRAN  PROGRAM  FOR  GENERATING  A 
_ mn 

The  ALTRAN  program  listed  in  this  appendix  is  used  for  computing 

the  polynomials  A  of  Eq.  (228).  A  recurrence  relation  similar  to 
mn 

Eq.  (153)  is  the  basis  for  the  routine. 

In  order  to  execute  this  program,  an  ALTRAN  compiler  is  required 
to  produce  MAP  and,  thereby,  the  final  object  code.f  The  accompanying 
routine  and  its  associated  control  statements  are  appropriate  for 
execution  on  The  RAND  Corporation's  IBM  7044  computer.  The  prologue 
to  the  listings  describes  the  program's  parameters,  usage,  and  limita¬ 
tions. 

t 

The  ALTRAN  compiler  consists  of  transfers  to  ALPAK,  a  group  of 
subroutines  for  operating  on  certain  polynomials  and  rational  functions. 
The  details  of  ALTRAN  and  ALPAK,  including  the  format  for  representing 
polynomials,  are  discussed  in  Ref.  38. 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

(CLOSE 
SIB JOB 
(FILE 
(FILE 
SItOIT 


THE  FOLLOWING  LISTING  IS  AN  ALTRAN  PROGRAM  FOR  COMPUTING 
THE  POLYNOMIALS  LAMBOAIMyN)  IN  THE  VAR  I ABLES  SIGMA  AND  BETA. 

THE  POLYNOMIALS  ARE  EXPRESSED  AS  CONSTANTS  TIMES  ISI&MA*MM-I  )  )X 
I8ETA«*I)  FOR  1*0, l». ..'N  ANO  FOR  M*l(2t...  ANO  FOR  N*l,2,...,H 
THE  LAMBOAIMyN)  ARE  ZERO  FOR  N  GREATER  THAN  M. 

THE  FOLLGWING  PROGRAM  CALCULATES  LAMBOAIM.N)  rOR  M*1 ,2*... . 10. 
THE  CONTROL  CAROS  LISTED  BELOW  INDICATE  THE  APPROPRIATE  DECK 
SET  UP  FOR  EXECUTION  ON  THE  IBM  704*.  BINARY  DECKS  ARE  NOT 
FULLY  LISTED.  THE  LAMBOAIMyN)  ARE  MADE  AVAILABLE  AS  BOTH 
PRINTEO  AND  PUNCHEO  CARO  OUTPUT. 


S.SU07, REWIND 
MAP, FILES 

*  S.FBIA* yNONr, *,BL0CK=10 
•S.FBOA* .NONE,*, BLOCK* 10 
U07.SRCH 


SIBLOR  TMG 
SI  BLOK  TMGIO 
SIBLOR  TMGDFN 
SIBLOR  ALTRAN 
SENTRY  TMG 

STORAGE  13000 

LAYOUT  !L'  SIGMA  18,  GAMMA  18 
ALGEBRAIC  (L)  LAMBDA! 10, 10) 

INTEGER  M, N, R , ! 

LAMBOAI 1 , 1 ) * l 

LAM80A(2,i:*2*SIGMA 

LAMBDA! 2 , 2 ) =-3*S IGNA+GAMMA 

PRINT  LAMBDA! 1,1) ,LAM8DA( 2,1) » LAMBDA! 2,2) 

PUNCH  LAMBOAI 1,1) .LAMBDA  12,1), LAM8DA! 2, 2) 

DO  20  M*3,  10 


I  *M- 1 


DO  15  N*l,I 

LAMBDA! M,N)=LAMBDA!M- 1, N)*!N*I SIGMA-GAMMA) ♦ I  *  I S IGMA+GAHMA) ) / ! M-N ) 


PRINT  M , N, LAMBDA ( M, N ) 

PUNCH  LAMBOAIM.N) 

15  CONTINUE 
lambda l M ,M ) = I 

00  16  R* l , I 

LAMBDA! M,M)=LAMBOAIM,M)*!M*!SI GMA-GAMMA )  ♦R*  ( SI GM A* GAMMA ) ) / l R-M) 

16  CONTINUE 

PRINT  LAMBDA ! M, M ) 

PUNCH  LAMBDAIM,M) 

20  CONTINUE 
STOP 
END 

FINISH 


SI3SYS 
(CLOSE 
SI  8 JOB 
SI  EDI T 


S.SU06, REWIND 
MAP, I00P1 
U06 , SRCH 


A1BMAP  AL  T RAN 

AIEOIT  U07,SRCH 

AIBLOR  ALFSRT 

AIBLOR  REAOF 

AIBLOR  RE ADO 

AIBLOR  REAOI 

AIBLOR  OUT 

AIBLOR  f -JfiCHP 

AIBLOR  ALF 

AIBLOR  ALP 

AIEOIT  IN 

AIBLOR  I  OREO  10/28/66 

ACOICT  IOREC 

C  BINARY  CAROS  DELETED 

A7EXT  IGRED 

C  BINARY  CARDS  OELETfcD 

ADKENO  I OREO 
AIBLOR  POSTXX  10/28/66 

ACOICT  POSTXX 

C  BINARY  CARDS  OELETfcD 

ATEXT  POSTXX 

C  BINARY  CARDS  DELETEO 

AOK.END  POSTXX 
AENTRY  ALT 

AIBSYS 

ACLOSt  S«SU07, REMOVE 
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Appendix  C 

A  FORTRAN  IV  PROGRAM  FOR  GENERATING  THE  COEFFICIENTS  IN  A 
_ _ mn 

The  FORTRAN  IV  program  listed  in  this  appendix  is  used  to  compute 

the  constants  in  the  polynomials  A  of  Eq.  (228)  .  The  routine  is 

mn 

based  on  Eqs.  (237),  (240),  and  (241)  of  Appendix  D. 

The  following  program  is  compatible  with  the  IBM  7044,  7094, 
and  360  series  FORTRAN  IV  compilers.  The  program  obviates  the  large 
storage  needed  in  the  ALTRAN  routine  discussed  in  Appendix  B.  The 
routine  also  allows  cross-checking  with  the  results  of  ALTRAN. 

The  prologue  to  the  listing  describes  the  parameters,  usage,  and 
limitations  of  the  program,  as  well  as  the  format  of  the  printed 
results.  All  of  the  computation  is  performed  in  double  precision. 


nr.  ooonooonoon 
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SIBFTC  COEFF 


THIS  ROUTINE  PRODUCES  THE  MATRIX  OF  LAMBDA  I  M,N)  COEFFICIENTS. 

M  RANGES  FROM  l  TU  15  ANO  N  RANGES  FROM  1  TO  M. 

THE  COEFFICIENTS  ARE  ZERO  FOR  ALL  N  GREATER  THAN  M. 

THE  COEFFICIENTS  OF  THE  VARIABLES  SIGMA**I  X  BETA**J  ARE  GIVEN 
FOR  EACH  POLYNOMIAL  LAMBDA  (  M ,  N )  ■  THE  OUTPUT  FORMAT  IS  AS  FOLLOWS 
M*  « N*  ,  POWER  OF  S IGMA=  ,  POWER  OF  BE TA=  ,  COEFFICIENT*  . 
ALL  COMPUTATION  IS  IN  OOU8LE  PRECISION.  FOR  M  GREATER  THAN  15 
HIGHER  PRECISION  IS  REQUIRED,  THE  INDEXING  MUST  BE  INCREASED 
AND  THE  DIMENSION  STATEMENTS  MUST  BE  ADJUSTED. 


DOUBLE  PRECISION  S ( 16, 1 6  I , F { 1 6 ) , AM , AN , A , ANS ,B *C 

INTEGER  R,P 

COMPLEX  Q 

LOGICAL  TEST 

WRITE(6,6) 

6  FORMAT! IH1////I0X,IHM,3X,  IHN, 3X , 1 AHPOWER  OF  SIGMA*3X, 

. 1 3HPOWER  OF  RETA,3X,30HC0EFFICIENT  OF  SIGMA-BETA  TERM  ////) 
DO  l  M*  l  ,  1 5 
S(M,M)=l.DO 
AM*DBL E(FLOAT(MJ  I 
SIM+l, 1 )=-AM*S(M,  l) 

Ml*M*l 

DO  1  K*M  1,15 
S ( M , K ) *0. DO 

1  CONTINUE 

DO  A  M*2,15 
AM*DBLE(FLOAT(M») 

DO  A  K*2 » 15 

A  S(M*1,KI=S(M,K-1)-AM*S(M,K) 

F ( 1 )=l.DO 
F ( 2 ) *F (  l  ) 

DO  2  I  *2  ,  l  A 

2  F(  I*l)*Fm*UBLE(FLOAT|  I)) 

DO  100  M=  1 ,  10 
AM*DBLE(FLOAT(M) ) 

M2*M-1 

DO  100  N -•*  l  ,  M 
AN*DBLEIFLOAT { N 1  I 
MN*M-N+ I 

A*( !-l.DO)**(M*NI 1 / ( F I N* 1 ) *F ( MN I J 
DO  100  11*1, M 
I S 1 GMA* I ABS ( M- I  I ) 

IF! I I.EU.M. ANO.N.NE.M)  GO  TO  100 

IBETA*IABS(M2-ISIGMAJ 

ANS*0.D0 

DO  3  K* l ,M 

MK*M-K  ♦  ! 

B-f-ANl •♦K*S(M,K J*F (MR) ♦FIK  5 
00  3  R* 1 ,MK 
MKR*M-K-R*2 
DO  3  P= 1 , K 


'■c*  *■*&*&■•*  WHMHmv**  1 
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KP-K-P+L 

IFIIM-R-PM).EG.  ISIGMA)  ANS*ANS*C /  (  f  ( R. )  *F  (  P  )  *F  {  MKR ) *F { KP ) ) 
3  CONTINUE 
ANS*-ANS*A 

0»CMPLX(0.,l.)**I8ETA 

1Q1*REAL(Q) 

IQ2*A!MAGIQ) 

TEST*! Q2.EQ.0. AND. IQl.LT.O 
IF!  TEST!  ANS—ANS 

IF ( TEST )  MRITEI6,7)  M,N, l SIGMA, IBETA, ANS 
7  FORMAT (9X,I2,2X,I2«9X, 12, 15X«I2«11X, 024*16) 

IF!  TEST  I  GO  TO  100 
TEST* IQ1*EQ*0* AND* l  Q2  *L  T. 0 
IF C  TEST )  ANS*-ANS 

IF  I  TEST )  HR  I TE ( 6,5)  M,N, ISIGMA, IBETA, ANS 
5  FORMAT  1 9X, I2,2X,I2,9X,I2, 15X, 12, ilX,D24. 16, 2X, IHI) 

IF f  TEST  I  GO  TO  100 

IF! 1Q1 .EQ*0. AND. IQ2.GE. 0)  WRITEI6.5)  M,N, ISIGMA, IBETA, ANS 
IF! IQ2.EQ.0.AND. IQl.GE.O)  WRITEI6.7)  M,N, ISIGMA, IBETA, ANS 
100  CONTINUE 
CALL  EXIT 
END 
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Appendix  D 

A  SERIES  REPRESENTATION  FOR  A 
_ mn 

For  reasons  discussed  in  Appendix  A,  a  series  representation  for 

A  is  a  useful  supplement  to  the  product  form  given  in  Eq.  (228) . 

Starting  with  Eqs.  (68^  and  (228)  for  a  and  A  ,  one  can  write 

mn  mn 
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Sub8titution  of  Eqs.  (234)  and  (235)  in  Eq.  (233)  therefore  results  in 


V  Iny  (-1)k4~P  smk)  "k("-fc)i(k-i)i 
Ann  n! (m-n)  1  /,/,/,  rip!  (m-k-r)  I  (k-p-1)  ! 

k-C  r*0  p=0 


.-'•"•‘(ifl)"'  (236) 


Since  “0,  the  series  representation  for  simplifies  to 


.  m  m-k+1  k 

in-^nfcor  l  <-1>1‘slk)  l  l 

k=  1  r=  1  p=  1 


-nP/nr+P 


rr  cr 


m-r-p+1  r+p-2 


( r- 1)  1  (p-1)  1  (k-p)  I  (m-k-r+1)  ! 


In  order  to  obtain  the  coefficients  c^  in  the  representation 


(237) 


,  V  tn-l-i  i 

Am  *  2.  ‘l  a  S 


(238) 


which  follows  from  Eq.  (237),  one  notes  that  r  and  p  must  be  selected 


according  to 


m-l-i  m-r-p+1 

o  ^  a 


l’1 


.  r+p-2 


r  =  1,2,.  .  , m-k+1 
P  =  1,2,...  ,k 
k  =  1 ,2 , . . .  ,m 


Thus,  for  any  value  of  i  and  r 


(239) 


m-l-i  =  m-r-p+1  i  =  0,l,...,m-l;  k  =  1,2,.. . ,m;  r  =  1 ,2 , . . . , m-k+1  (240) 

so  that  p  is  constrained  to  the  value(s) 


p  -  i-r+2  with  p  c  (l,2,...,k) 


(241) 
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Consequently ,  is  the  aggregate  of  all  terms  of  Eq.  (237)  obtained 
by  allowing  k  to  range  over  the  integers  from  1  to  m;  r,  from  1  to 
m-k+1;  and  p  =  i-r+2,  with  p  contained  in  the  set  k  =  Such 

a  summation  is  easily  programmed  for  digital  computation;  a  FORTRAN  IV 
routine  for  finding  the  c^  is  provided  in  Appendix  C. 
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Appendix  E 

A  FORTRAN  IV  PROGRAM  FOR  COMPUTING  EXPONENTIAL  APPROXIMATIONS  g(t) 

The  accompanying  program  is  designed  to  calculate  the  exponential 

A  2 
approximation  g(t)  to  a  prescribed  function  g(t)  c  L  (0,®).  The 

algorithm  is  based  on  the  recurrence  relations  Eqs.  (157)  and  (158) 

for  the  m^  orthonormal  basis  function  X  ( t)  ,  on  recurrence  relation 

m 

Eq.  (153)  for  X  ,  and  on  Eqs.  (132)  or  (162)  for  the  moments  of  g(t) . 
mn 

These  moments  are  approximated  by  a  64-point  Gaussian  quadrature 
scheme  and  are  used  to  form  the  Fourier  expansion  coefficients  a 

m 

A 

according  to  Eq,  (131).  The  approximant  g(t)  is  finally  obtained 
from  Eq.  (133)  . 

The  numerical  examples  discussed  in  Section  IX  have  been  solved 
with  the  use  of  program  APRX.  The  values  of  o  and  3  obtained  from 
Eqs.  (175),  or  from  Eqs.  (184)  and  (185),  are  input  parameters  to 
the  program  as  described  in  the  prologue  of  APRX.  All  other  program 
options,  definitions,  and  limitations  are  also  clarified  in  the 
listings . 

The  real  and  imaginary  parts  of  g(t)  must  be  supplied  as  the 
double-prr ’ ision  function  subroutines  named  GR(T)  and  JI(T) ,  respec¬ 
tively,  Once  APRX  is  entered,  GR,  G1 ,  and  all  the  supporting  routines 
are  automatically  invoked  to  produce  a  51-point  tabular  and  graphical 

A  A 

display  of  g(t).  Card  output  of  g(t)  for  subsequent  pro'essing  is 


also  available. 


onooonc-inooooononooooooooooooooooor 
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tlBFTC  APRX 
C 


THIS  ROUTINE  IS  THE  MAIN  PROGRAM  FOR  OBTAINING  AN  EXPANSION  OF 
A  GIVEN  FUNCTION  GIT)  IN  TERMS  OF  THE  ORTHONORMAL  SET  X ( M*  T  )  . 

THE  REAL  PART  OF  GIT)  MUST  BE  SUPPLIED  AS  A  DOUBLE  PRECISION 
FUNCTION  SUBROUTINE  CALLED  GRIT).  AN  EXAMPLE  IS  GIVEN  AT  THE 
ENO  OF  THE  LISTINGS.  THE  IMAGINARY  PART  OF  GIT)  MUST  BE  SUPPLIED 
AS  A  DOUBLE  PREC1SON  FUNCTION  SUBROUTINE  CALLED  GIIT).  AN  EXAMPLE 
IS  GIVEN  AT  THE  END  OF  THE  LISTINGS.  THE  REPRESENTATION  IS  BEST 
IN  AN  INTEGRAL  SQUARE  ERROR  SENSE  OVER  THE  INTERVAL  ZERO  TO 
INFINITY.  IN  ORDER  TO  FACILITATE  THE  INTEGRATION  AND  PLOTTING 
THE  FOLLOWING  PARAMETERS  MUST  BE  READ  IN  FROM  A  DATA  CARD 
ACCORDING  TO  THE  FORMAT  NUMBER  8  GIVEN  BELOW. 

MC*THE  NUMBER  OF  TERMS  DESIRED  IN  THE  EXPANSION. 

SIG-THE  VALUE  CHOSEN  FOR  SIGMA. 

BET-THE  VALUE  SELECTED  FOR  BETA. 

AA«TH£  SMALLEST  VALUE  OF  T  I  GREATER  THAN  OR  EQUAL  TO  ZERO)  BEFORE 
WHICH  THE  PRESCRIBED  FUNCTION  GIT)  IS  ESSENTIALLY  ZERO. 

BB-THE  LARGEST  VALUE  OF  T  I  GREATER  THAN  ZERO)  AFTER  WHICH  GIT) 

IS  ESSENTIALLY  ZERO.  AA  AND  BB  ARE  THE  LIMITS  USED  IN  COMPUTING 
THE  FOURIER  COEFFICIENTS  FOR  THE  EXPANSION  ON  ZERO  TO  INFINITY. 
GRHAX«MAXIMUM  VALUE  OF  THE  REAL  PART  OF  GIT)  TO  BE  IN  THE  GRAPH. 
GRMIN«MINIMUM  VALUE  OF  THE  REAL  PART  OF  GIT)  TO  BE  IN  THE  PlOT. 

THE  SAME  ORDINATE  SCALE  WILL  BE  USED  IN  PLOTTING  THE  IMAGINARY  PART. 
TMAX*MAX IMUM  ABSCISSA  VALUE  TO  BE  USED  IN  THE  PLOT  OF  GIT). 
TMIN^MINIMUM  ABSCISSA  VALUE  TO  BE  USED  IN  THE  PLOT  OF  GIT). 

IPUNCH-0  IF  THE  PLOTTEO  POINTS  OF  GRIT)  AND  GIIT)  AND  THE 
APPROXIMATE  VALUES  OF  GHATRIT)  AND  GHAT  I ( T )  ARE  NOT  TO  BE  PUNCHED 
ON  DATA  CARDS  (ACCORDING  TO  FORMAT  NUMBER  1  BELOW). 

ALL  COMPUTATION  IS  IN  DOUBLE  PRECISON.  !F  MC.  THE  NUMBER  OF 
TERMS  IN  THE  EXPANSION,  IS  GREATER  THAN  20,  THE  PROGRAM  DIMENSION 
STATEMENTS  MUST  BE  ADJUSTED.  ALL  SUBSEQUENT  SUBROUTINES  ARE 
AUTOMATICALLY  CALLEO  ONCE  THE  MAIN  ROUTINE  APRX  IS  ENTERED. 


DOUBLE  PRECISION  CFRI 20 ) ,CF 1 1  20 ) ,GHATR, GHAT  I , T * S IG , BET 
. ,ERROKR,KMSE,AA,B8,S2P,ESIGT,Pl , 
•LSR(20,20),LSI(20,20),S1(20),S2(20) 

COMMON  /COEFF/  CFR,CFI 
COMMON  /MOMENT /SI « S2 
COMMON  /LSCOM/LSR.^S! 

COMMON  S IG, BE  T , MC , NNN, AA, BB , S2P,ESIGT,T  »GHATR,GHAT I 
EXTERNAL  ERRORR 

REAL  GRAPH! 1000 ) , OROR 1 5 1) ,OROI ( 5 1 1 , ABSC ( 5 1 1 ,H( 51 1 ,HH< 5 1 1 
DATA  PI/).  14159265358979300/ 

S2P*DSQRT! 2.D0*PI 5 

l  RE AD  I  5 , 8  I  MC,SIG»PET*AA,BB,GRMAX,GRMIN, TMAX , TM IN, I  PUNCH 
8  FORMAT! I ),2D16.  1,  205. I , 4F 5. I , I  5  ) 

WRITE(6,23)  S1G,BET,MC,AA,BR 

23  FORMAT! 1 H 1 , 6HS I GMA* , 02*. 16//6H  BE *A«,D24. 16//8H  FOR  M-  ,I2,6H  TERM 
.$  //9H  FROM  T *  , D24. 1 6/ / 7H  TO  T»  ,024.16//) 

CALL  LAMBDA 
WRITE  16, 18) 

18  FQRMAT1//27H  THE  MATRIX  OF  LAMBDA! M,N)  //) 
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DO  19  M* 1 »  MC 

19  WRITEI6.20)  I M,N, LSR ( M, N) ,LS I I M( N) ,N« I ♦ M ) 

20  F0RMATI3H  LI ,  1 2  ♦  IH,  12  ,  >H)  * ,  024.  1 6, AH  *1  ,024.16) 

WRITE (6 , 17 ) 

17  FORMAT I // 1 0H  THE  MOMENT  VECTOR  //) 

CALL  COEF 

WRI tC (6, 16)  <I,Sl<I),S2II),I*l,MC) 

16  FORMAT  1 3H  G(,I2»2H)*,D24.16,4H  *1  ,D24.16  I 
WR I TE<6 « 900) 

900  FORMAT { ///70H  THE  REAL  AND  IMAGINARY  PARTS  OF  THE  EXPANSION  COEFF 

•  IC I  ENTS  ,  AIM)  ///) 

WRITE(fc»901 )  I  I,CFR(Ji)»CFl(  I  ) « I  *  l,MC ) 

901  F0RMATI4H  A( ,  I2,2H)*,D24.16,5H  +  J  ,024.16) 

DlV*ITMAX-TMIN)/50. 

DO  10  1*1*51 

ABSC I  1 1 *AA*FLOAT  t l-l ) *01 V 
T«DBLE(^8SC( I ) ) 

HI  1 )*GR(T) 

HH<  1 )  *G  1  (  T ) 

CALL  GHAT 
ORDRI 1)*GMATR 
ORDI II) *GMAT I 
I c I  I  PUNCH )  11,10,11 

11  WRITEI7,7)  ABSCI I), HI i ) ,HH( I ), OR OR  I  I ) *0RDI I  I ) 

7  FORMA! I 5F 10. 4 ) 

10  CONTINUE 

CALL  INTGRL (AA,BB»ERRORR,RKSE) 

RMSER*OSQRT I RMSE ) / I BB-AA) 

CALL  PL0T2!GRAPH,TMAX,TMIN,GRMAX,GRM!N) 

CALL  PL0T3I 1H>, ABSCI 1I»H( 1) ,51) 

CALL  PL0T3I1H*, ABSCI 1 ) , ORDR I l ) , 5 l ) 

WR I TE 1 6  *  008 ) 

808  FORMATI /// ) 

WR I TE 1 6 , 800 )  I ABSCI I i , ORDRl I ) ,ORDI I  I ),H(I),HHl I ) *1*1,51  ) 

800  FORMAT  1 4H  T* , F6. 2 e 3X , 9HRE I  GHAT ) *, E 17. 8, 3X, 9H I  MI  GHAT ) *, E 1 7. 8 , 3X , 

. 6HRE (G) *, E 17.8, 3X,6HIM1 G)*,E 17.8  ) 

WRITE  16, 812)  RNSER 

812  FORMAT I////15H  THE  RMS  ERROR*  ,  E10.3/) 

WR I TE I  6 , 809  ) 

809  FORMATI  1HL) 

CALL  PL0T4I l , IH  ) 

WR I TEf 6,801 ) 

801  FORMATI 35H0**RE I GI T ) )  **RE  t  GHA  TIT))  ) 

CALL  PL 0T2I GRAPH, TMAX, f MI N, GRMAX , GRM! N) 

CALL  PL0T31 1H+, AB'X I  1 ) , HHI 1 ) ,5 1 ) 

CALL  PL0T3I IH*, ABSC I l ) , ORDI I l ) ,5 1 5 
WRITEI6,809) 

CALL  PL0T4I 1 , IH  ) 

WR I TE 16,802 ) 

802  FORMATI 3SM0** IMI Gl T ) )  *■ I M (GHAT IT))  ) 

GO  TO  1 

CALL  EXIT 
END 


ooooonow*  o  o  <  >  o  rv  o  o  *► 
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IBFTC  LAMBDA 


THIS  SUBROUTINE  COMPUTES  THE  CCEFF ICENT S  LAH80A(M,N)  IN  DOUBLE 
PRECISON.  THE  REAL  AND  IMAGINARY  PARTS  OF  LAMBOA(M,N)  ARE  STORED 
IN  ARRAYS  LSRIM.Nt  AND  LSIIM.N)  RESPECTIVELY. 


SUBROUTINE  LAMBDA 

DOUBLE  PRECISION  ISR(  20, 20 ) ,L S l( 20 « 20) , S I G,8E i , PI • A,B,C, AM, AN, 
. Q«  E  ,  F  ,G,H,AK,D 
OATA  PI/3.1M592653589793/ 

COMMON  SIC, BET, MC 
COMMON  /LSCOM/LSR.LSl 
A»SIG*SIG4BET*BET 
B«SIG*$iG~BET*BET 
C*2.D0^SIG*BET 
LSRII»1)-0SQRT(S!G/PI ! 

LSI ( l, l 1*0.00 
DG  l  M= I , NC 
AM  =  M 

E*(-l)**H 

F*C*AH 

G*B*AM 

0*CSCRT ( ( AM+l .00 ) /AM ) 

00  ?  N= 1 , M 
AN*N 

c*o/(a»(mh+i.lg-an) : 

LSR ( M* 1 , N ) * ( l SR(M,N)*( AN*  A*  G)-LSl|H,N)*F)*8 

2  LSI  lM*lfN|*(L$I(M,N)*IAN*A+  G)*-LSR(N,N)*F)*Q 
H-OSQRTI  AMM.DO) 

F*E*H*LSR| 1,1) 

G*0. DO 

00  3  K*  l  *  M 

F*F-LSR(M*1,R) 

3  G«G-LSI(M*l,K) 

LSR(M*l,M*l)=F 

1  LS 1 1 M* 1 » M+ 1 ) =G 
RETURN 
END 


IBFTC  COEF 


THIS  SUBROUTINE  COMPUTES  THE  FOURIER  COEFFICIENTS  AIM)  IN  DOUBLE 
PRECISON.  THE  REAL  ANO  IMAGINARY  PARTS  OF  THE  RESULTS  ARE 
STOREO  IN  THE  ARRAYS  CFR(M)  ANO  CFHM;  RESPECTIVELY. 


SUBROUTINE  COEF 
tXTERNAL  RINTG, I INTG 

DOUBLE  PRECISION  I  I NTG,  R I  NT G, S I G , BET , S2P , LSR 1 20, 201 , AA, BB 


oonnnoono** 
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• *LS I (20,20) .CFRC201.CFI (20) . SI 1 20! . S2I 201 .P 1 ,C 1 1  20,20) ,C2I 20,201 
..C3(20,2CI,C*‘20,20I .ESIGT.T 
t-OMHON  /COEFF/CFR.Cf I 
COMMON  /LSCOM/LSR.LSI 
COMMON  /ARRAY /CI.C2.C3.CA 
COMMON  /M0MENT/S1.S2 

COMMON  S1G,8£T,MC.NNN,AA,BB,S2P.ESIGT,T,GHATR,GFATI 

00  10  N* i «  MC 

NNN=N 

CALL  IM  jKl  l  AA  ,  ;>R  »  R  l  NTG  ,  S  )  INI) 

CALL  INTuHL(AA,UH, I INTG.S/IN) ) 

1C  CONTINUt 
00  1  M*  l , MC 
CFR  C  Ml =O-D0 
CFI (M1*0.00 
CO  S  N* 1 »M 

CFR(M|*CFRCM)*LSRCM,N1*S1 tN)*LSl (M,NI*S2(N) 

CFI (Ml =CF I(MJ*LSR(M,N)*S2(N)-LSI (M,N)*S1(N) 

5  CCNTINUE 

CFRCMJ*S2P*Cf RIM) 

CFHM)  =  s2P*rr;tM) 

1  CONTINUE 

00  ?  M=1,MC 
DO  2  N=  1  » M 

C1«N,N)=CFR(M)*LSP(M.N)-CFI (M)*LSIIM,N) 
C2IMtNl=-CFR(M)«LSHM.Nl-CFHM)*LSR(M,N) 

L3IM,N1*CFI CM)*LSRIM,N}*CFRCMJ*LSICM,N) 

2  CAIM.N)=-CFI  CM}*LSUM.N1*CFR(M|*LSRCM,NI 
RETURN 

END 


1BFTC  GHAT 


THIS  SUBROUTINE  FORMS  THr  REAL  ANO  IMAGINARY  PARTS  OF  THE 
APPROXIMATION  GHAUT)  TO  THE  PRESCRIBED  FUNCTION  GIT).  ALL  THE 
COMPUTATION  IS  IN  DOUBLE  PRECISION.  THE  REAL  AND  IMAGINARY  PARTS 
OF  THE  APPROXIMANT  ARE  STORfcD  IN  GHATR  AND  GHAT  I  RESPECTIVELY 
FOR  EACH  VALUE  OF  THE  INDEPENDENT  VARIABLE  T. 


SUBROUTINE  GHAT 

DOUBLE  PRECISION  C l< 20, 20 J , C2 i 20, 20 ) ,C 3 ( 20, 20 J ,C M 20 , 20) , EE ,0E 
.  ,GHATR,GHATI,T,S IG, BET, AAt BB.S2P.ES IGT,eT ,CC,SS 
COMMON  /ARRAY/C 1.C2.C3.CA 

COMMON  SIG.BET, MC ,NNN, AA . BB.S2P.ESIGT.T ,GHA I R .GHAT  I 
BT«BET*T 

ESIGT«DEXP{-SIGYT) 

GHATR»O.DO 
GHAT  1*0.00 
00  l  M=l,MC 
OE* l .00 
00  l  N* l  ,H 


o  o  o  o  o  nonoooo*  oooooon 
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C£*DE*ESIGT 
£E»8T*DBLE  I  FLOAT  I N) ) 

CC*DC0S ( EE ) 

SS-DSIN(EE) 

GHATR«GHATR*DE*(CC*CllMtN)*SS*C2(MtN)  ) 
I  GHAT  I -GHAT  I *DE* ( CC*C J< M , N) ♦ S$*C4 ( M , N ) ) 
GHATR*GHATR*S2P 
GHAT 1*GHAT I  *S2P 
RETURN 
END 


AIBFTC  RINTG 


THIS  SUBROUTINE  COMPUTES  THE  REAL  PART  OF  THE  MOMENTS  OF  THE 
PRESCRIBED  FUNCTION  GIT).  THESE  ARE  USED  IN  COMPUTING  THE 
FOURIER  EXPANSION  COEFFICIENTS.  COMPUTATION  IS  IN  DOUBLE  PRECISON. 


DOUBLE  PRECISION  FUNCTION  r INTGf  T) 

DOUBLE  PRECISION  GR.GI , S I G, BET , AN, T 

COMMON  SIGtPET tMC«N 

AN*N 

RINTG*DEXPl— AN*S IG*T  )  *!  GR I  T  )*DCOSI AN*BET*T  )«-GIIT)*DSlNl  AN*BET*T  )  ) 

RETURN 

END 


1BFTC  IINTG 


THIS  ROUTINE  COMPUTES  THE  IMAGINARY  PART  OF  THE  MOMENTS  OF  THE 
PRESCRIBED  FUNCTION  GIT).  THESE  ARE  USED  IN  COMPUTING  THE 
FOURIER  EXPANSION  COEFFICIENTS.  COMPUTATION  IS  IN  DOUBLE  PRECISON. 


DOUBLE  PRECISION  FUNCTION  UNTGIT) 

DOUBLE  PRECISION  GR,GI, SIG,BET, ANtT 
COMMON  SIG* RET  »MC  »N 

AN»N 

IINTG»DEXPI-AN*SIG*T)*I-GRI T ) *DS INI AN*BET*T)*GI !T)*DCOSI AN*BET*U) 

RETURN 

END 


IBFTC  INTGRL 


THIS  PROGRAM  PERFORMS  INTEGRATION  IN  DOUBLE  PRECISION  AND 
IS  BASED  ON  THE  6*  POINT  GAUSSIAN  QUADRATURE  FORMULA. 
A«LOWER  LIMIT  OF  INTEGRATION  IN  DOUBLE  PRECISION. 


o  n  n  o  o 
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B-UPPER  LIMIT  OF  INTEGRATION  IN  DOUBLE  PRECISION. 

F*THE  NAME  OF  THE  OOUBLE  PRECISION  FUNCTION  TO  BE  INTEGRATED. 
ANS*THE  RESULTANT  INTEGRATION  IN  DOUBLE  PRECISION. 


SUBROUTINE  I NTGRL ( A, B ,F , ANS I 
EXTERNAL  F 

OOUBLE  PRECISION  A,B,F, ANS* X( 32  I ,H( 32) t SI » S2 fUl28l • VI 4) , Y( 28 ) , 2 ( A) 
DATA  U/.024350292663424432509, .072993121787799039450,. 121462819296 
11205544  70,.  1696444204239928 1803 7,. 2 1 742 3643 740007084 150,. 264687 1 62 

2208767416374..  31 13228 71  *902 10956 1  58 ,. 35 72201 58 3376681 15950, .4022 70 

31 5 7963991603696..  446366 01 7253464087985. .489403 14570 7052957479. . 5 3 1 

4279464019894545658. . 571 89 564620 2634 034284. .6111553551 72393250249 , • 

564896547 1254657 339858..  6852 363 13054233242564.. 71988 1850 1 7 161082684 

69..  7528 1990726053 189661 2.. 78397235894 334 140 76 10.. 81 32653 15 12279755 

79742..  840629296252580362752.. 8659993981 5409281 9761. .88931544599511 

84 105 85 3..  9 105221 37078502805756. .929569172131939575821.. 9464 1137485 

984028 16062..  961008799652053 7 1 89 19.. 9 7332682 7789910963742/ 

DATA  V/  .98333625 

103884625956931..  99 101 33 7 1476744320739.. 9963401 16 77 19552 79347.. 9993 
1105041735772139457/ 

OATA  Y/. 0486909570091 39720383,. 0485 7546 7441 503426935,. 048344762234 

18029571 70..  047999388596458307728.. 0475401657 14830308662.. 0469681 82 

281621 00 17325. . 04628479b 58131441 7296, . 0454916279274181 44480, .044590 
3558163756563060,-043583724529323453377, .042473515123653589007, .041 

42625632 42 62 35 2 86 10. . 039 9537411 32 720 34 1387. .03855015317861562 9129.. 

5037055 128540240046040..  0354722 132568823838 11 ..03380516183714160939 

62..  032057928354851 553585.. 030234657072402478868.. 02833967261425948 

73228..  02637746971 5054658672.. 024352702568710873338.. 0222 701 7380838 

83254159..  0201 34823153530209372.. 01 795 17 15775697343085.. 01572603047 

960247 19 322..  01 34630478967 18642598.. 01 1 168 1 394601 3 1 1288 19/ 

DATA  It  .00884675 

10982636394772  3, .006504 4  57968978362856,. 00414  70  33260562467635, .0017 
1183280721696432947/ 

OATA  I  START/* 1 / 

Si* ( B-A 1 /2.D0 
S2*(B*AI/2.00 
I F ( ISTARTI  4,4,5 
5  DO  2  1*1,28 
X{ I 1*U( I  I 

2  W< I  1 *Y I  I ) 

DO  3  1*1,4 
X( I *28 1  —  V <  I  1 

3  W( I *28 1 *Z (  I  ) 

4  ANS*0. DO 

00  1  1*1,32 

1  ANS*ANS*W{ I)*(FIS1*X(I)*S2)*F(-S1*X(I)*S2I) 

ANS*ANS*Sl 
1  ST  ART  *-0 
RETURN 
ENO 


IIBFTC  ERRORR 
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C 

c 

C  THIS  ROUTINE  COMPUTES  THE  INTEGRAL  SQUARE  ERROR  BETWEEN  THE 

C  PRESCRIBED  FUNCTION  ANO  THE  APPROXIMANT  GENERATED  BY  THE 
C  ACCOMPANYING  PROGRAMS.  COMPUTATION  IS  IN  DOUBLE  PRECISION. 

C 

C 

DOUBLE  PRECISION  FUNCTION  ERRORR ( X) 

DOUBLE  PRECISION  S IG,BE T, AA, BB« ESIGT, S2P, T , X ,GHATR,GHAT I ,GR,G I 

COMMON  S IG, BET ,  MC«  NNN, AA, BB,  S2P»ESIGT»T ,  GHATR  ,-GHAT I 

T*X 

CALL  GHAT 

ERRORR* I GHATR-GK I  X ) ) **2 ♦ I  GHAT I-G  I  ( X  » ) **2 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  XR(M.T) 

HBFTC  XR 

C 

C 

C  THIS  ROUTINE  COMPUTES  THE  REAL  PART  OF  THE  M  TH  ORTHONORMAL 
C  BASIS  FUNCTION  X(M,T).  COMPUTATION  IS  IN  DOUBLE  PRECISION. 

C  THIS  PROGRAM  IS  USED  ONLY  FOR  CHECKING  ORTHONORMALITY  AND  IS 

C  NOT  OTHERWISE  USED  BY  THE  ACCOMPANYING  PROGRAMS. 

C 

C 

DOUBLE  PRECISION  LSR120.20) ,LSII20,20> ,S2P.AN,SIG,BET, T,PI 
. .ESIGT, OE, EE, AA,BB 
COMMON  /LSCOM/LSR.LSI 
COMMON  SIG*BET , MC * NNN » A A , BB , S2P, ESIGT 
DATA  PI/3. 141592653589703/ 

ESIGT*DEXP(-SIG*T) 

DEM.  DO 
XR-0.00 
DC  l  NM,M 
AN*N 

DE*DE*E  S IGT 
EE*AN*8ETPT 

1  XR«XR*OE*I LSR I M#  N) *DCOS (EE)-LSI(MfN)*DSIN(EE)) 

XR*S2P*XR 

RETURN 

END 


$1 BFTC  XI 
C 

c 

c  THIS  ROUTINE  COMPUTES  THE  IMAGINARY  PART  OF  THE  M  TH  ORTHONORMAL 

C  BASIS  FUNCTtON  X(M,T).  COMPUTATION  IS  IN  DOUBLE  PRECISION. 

C  THIS  PROGRAM  IS  USED  ONLY  FOR  CHECKING  ORTHONORMALITY  AND  IS 

C  NOT  OTHERWISE  USED  BY  7HE  ACCOMPANYING  PROGRAMS. 

C 
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DOUBLE  PRECISION  FUNCTION  XllM,T) 

DOUBLE  PRECISION  L SR ( 20, 20 ) , L S II 20  20 ) , S2P , AN, S IG, BET , T , P  I 
., ESIGT.DE, EE, AA,8B 
COMMON  /LSCOM/LSR.LSi 

COMMON  SIG,BET,MC,NNN,AA,BB,S2P,ESIGT 
ESIGT*DEXP(-SIG*T) 

DE*1.D0 
X I *0. DO 
DO  l  N* 1 ,M 
AN*N 

£E»AN*BET*T 

DE«DE*E$IGT 

l  XI*XI*OE*lLSR(M»N)*DSINIE£)+LSI ( M,N ) *DC OS l EE ) ) 

XI«S2P*XI 

RETURN 

END 


IIBFTC  GR 


THIS  ROUTINE  SUPPLIES  THE  REAL  PART  OF  A  SAMPLE  FUNCTION  GIT). 
T»THE  INDEPENDENT  VARIABLE  IN  DOUBLE  PRECISION. 

GR*THE  REAL  PART  OF  GIT)  IN  DOUBLE  PRECISION. 


DOUBLE  PRECISION  FUNCTION  GRIT) 
DOUBLE  PRECISION  T 
I F I T . LE . 0. DO )  GO  TO  I 
IFIT.LE..5D0)  GO  TO  2 
IFIT.LE..95D0)  GO  TO  3 
GR»0EXP(-2.*2377D0*T) 

GO  TO  10 

1  GR»O.DO 
GO  TO  10 

2  GR-2.D0*T 
GO  TO  10 

3  GR»-2.D0*IT-l.D0) 

10  RETURN 

END 


IBFTC  GI 


THIS  ROUTINE  SUPPLIES  THE  IMAGINARGY  PA»r  OF  A  SAMPLE  FUNCTION  GIT). 
T»THE  INDEPENDENT  VARIABLE  IN  DOUBLE  PRECISION. 

GI-THE  IMAGINARY  PART  OF  GIT)  IN  DOUBLE  PRECISION. 


DOUBLE  PRECISION  FUNCTION  GI(T) 
DOUBLE  PRECISION  T 
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Gl»0.00 

RETURN 

END 


IIBNAP  PLOT 

THIS  PLOTTING  ROUTINE  IS  AN  ADAPTATION  OF  SHARE  UHPLOT  FOR  THE 
IBM7044  .  THE  ROUTINE  UTILIZES  THE  FORTRAN-MAP  INPUT-OUTPUT 
r 'OGRAM  FACILITY  WHICH  IS  SUPPLIED  AS  A  SEPARATE  SUBROUTINE 

♦  WRITES  GRAPH  IMAGES  ON  OUTPUT  UNIT  G 

*  CALLING  SEQUENCES  ARE 


ENTRY 

ENTRY 

ENTRY 

ENTRY 

ENTRY 

ENTRY 

ENTRY 


CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

PLOT  I 

PLOT  2 

PLOT3 

PL0T4 

FPLOT  4 

OMIT 

PLTAPE 


PLOT  1 ( NS CAL E  » NH*  SBH*NV , SBV ) 
PLOT2I IMAGE, XMAX,XM IN, YMAX,YMIN) 
PLOT3(8CD»X,Y «NDATA ) 

PLQT4 ( NCHAR  *NHABCDEF ...  ) 
OMITIARG) 

PLTAPE 1 1  TAPE  I 


EXTERN  WDATA 


* 

* 

COL 

SPACS 


* 

PLOT  1 
* 


COLUMNS  IN  OUTPUT  LINE  (14011 
SET  UNUSED  SPACES  AT  RIGHT  EDGE  OF  PAGE 
OR  CARD.  SPACS  MUST  BE  AT  LEAST  6 


EQU  COL-ll-SPACS  COLUMNS  IN  OUTPUT  LINE  AVAILABLE 
SPACE  5 
PLOT  l 

REM  MAIN  JOB  OF  PL0T1  IS  TO  EXAMINE  ARGUMENTS  AND  PREPARE 

REM  SAMPLE  GRIDLINE  I  DASH  TO  DASH-WORDSU )  AND  SAMPLE 

REM  NON-GRIO  LINE  (BLANK  TO  BL ANK-WOROS+l )  FOR  PL0T2 


IMAGE 


SAVE 

ClA 
STA  [ 
AOO 
STA  ( 
ST  Z  k 
CLA  C 
STO  k 
CIA* 


ENTRY 


PLOTl 


3.4 
DELTA 

FIVE 

DELT 

WRONl 

ONE 

WR0N3 

4.4 


SCALE  FACTORS  AND  DECIMAL  POINT  POSITIONS 


WRONl 


CLEAR  ERROR  FLAG,  PLOTl 


TSX  FIX, 2 
TZE  ERR l 
STO  NH 
CLA*  5,4 
TSX  FIX, 2 
TZE  ERR l 


WR0N3  *  1  SET  ERROR  FLAG,  MISSING  PLQT2 

NH,  NUMBER  OF  HORIZONTAL  GRID  LINES 

ZERO  ARGUMENT  ILLEGAL •  ERROR  RETURN 

SBH,  NO.  OF  SPACES  BETWEEN  HORIZ.  GRID  LINES 

ZERO  ARGUMENT  ILLEGAL,  ERROR  RETURN 
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STO  SBH 
IDQ  SBH 
MPY  NH 
STQ  LINES 
CIA*  6,4 

TSX  FIX, 2 
TZE  ERR  1 
STO  NV 
CLA*  7,4 

TSX  FIX, 2 
TZE  ERRI 
STO  SBV 
LOO  SBV 
MPY  NV 
STO  TOT 
CLA  TOT 

AOO  ONE 
STO  TOTAL 
SUB  CHID 
TMI  PASS 


LINES  *  NH*SBH  MAXIMUM  LINE  INDEX 

NV,  NUMBER  OF  VERTICAL  GRID  LINES 

ZERO  ARGUMENT  ILLEGAL,  ERROR  RETURN 

SBV,  NO.  OF  SPACES  BETWEEN  VERT.  GRID  LINES 

ZERO  ARGUMENT  ILLEGAL*  ERROR  RETURN 


TOT  «  $BV*NV 


MAXIMUM  COLUMN  INDEX 


TOTAL  *  TOT  ♦  I  TOTAL  COLUMNS  PER  LINE 

WHENEVER  TOTAL  .G.  GRAPH  WIDTH,  ERROR  RETURN 


* 

ERRI 


CAL 

CALL 

PZE 

PZE 

PZE 

PZE 

CLA 

STO 

RETURN 


OTAPE 
WOATA 
FORM 
ER I «  0 , 1 
WRONG, 0, 
0 

FPONE 
WRON1 
PLOT  1 


RETURN  L,  UNSUCCESSFUL  PLOTl 
WRON1-1,  SET  ERROR  FLAG,  PLOTl 
UNSUCCESSFUL  PLOTl 


PASS  CLA  TOTAL 

TSX  FLOAT, 2 
FOP  SIXF 
STQ  TEMP 
CLA  TEMP 
FAO  N999 
TSX  F  I X I , 2 
STO  WORDS 
LOP  WORDS 
MPY  SIX 
STQ  TOTLS 
LXA  WORDS, 2 
CLA  TOTLS 
SUB  TOTAL 
PAX  0,1 
CLA  OSH, 1 
STO  0ASH«1,2 
LOQ  BLNKK 
STQ  BLANK* 1,2 
TIX  GA1 ,2, 1 

TRA  GA2 


WORDS  =  TOTAL/6.  ROUNDED  UP  TO  NEAREST  INTEGER 
WORDS,  NUMBER  OF  MACHINE  LOCATIONS  PER  LINE 


TOTLS  *  W0R0S*6 


BCD  CHARACTERS  PER  LINE 


LAST  WORD  OF  A  HORIZONTAL  GRID  LINE 

SET  UP  LAST  WORO  IN  HORIZONTAL  GRID  LINE  IMAGE 

LAST  WORD  OF  NON  GRID  LINE 

SET  UP  LAST  WORD  IN  NON  GRID  LINE  IMAGE 

ONE  WORD  PER  LINE  CASE 


o  o  o  o  o  o  o  o  o  o  o 
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GAI 

CLA 

OSH 

G 

GA3 

STO 

DASH* 1,2 

SET  REMAINDER  OF  HORIZ.  GRID 

G 

STQ 

BL ANK> 1 • 2 

SET  UP  REMAINDER  OF  NON  GRID  LINE  IMAGE 

T  I  X 

GA3,2, 1 

G 

• 

G 

GA2 

ST  Z 

I 

COL.  INDEX  FOR  VERTICAL  GRID 

G 

* 

G 

GAMMA 

T$x 

PLACE, A 

PUT  VERTICAL  GRID  1  IN  IMAGE 

G 

PZE 

I  EYE 

PZE 

l 

PiE 

blank 

* 

G 

TSX 

PLAC8,A 

PUT  ♦  AT  INTERSECTIONS 

G 

PZE 

1  PLUS 

PZE 

1 

PZE 

DASH 

* 

G 

CLA 

I 

ADD 

SBV 

STO 

I 

I-I+SBV,  INCREMENT  COLUMN  INDEX  FOR  VERT  GRID 

SUB 

TOTAL 

TZE 

GAMMA 

IF  ZERO  OR  MINUS  LINE  IS  UNFINISHED,  RETURN 

TMI 

GAMMA 

* 

G 

DELTA 

CLA 

** 

NSCALE 

NSCALE,  DETERMINES  SCALE  FACTOR  MODIFICATION 

TZE 

ETA 

STANDARD  SCALE  FACTORS  AND  DEC  POINT  POSITIONS 

* 

G 

AXT 

A»  A 

GT*  OEC  POINT  POSITION  FOR  X 

DELT 

CLA 

♦  ♦,A  NSCALE+5  G5  *  SCALE  FACTOR  FOR  X 

G 

TSX 

TF I X  t  2 

GA*  OEC  POINT  POSITION  FOR  Y 

STO 

G3+A  t A 

G3  *  SCALE  FACTOR  FOR  Y 

G 

TIX 

CELT , A* 1 

* 

G 

ETA 

CLA 

GA 

TZE 

GAT 

G 

TPL 

GAT 

G 

ZAC 

NEG.DEC.PT.  POSITION  *  0 

G 

TRA 

GAS 

G 

GAT 

CAS 

EIGHT 

MAX.  DEC.  PLACES,  ORDINATE,  *8 

G 

CLA 

EIGHT 

IF  GA  GTK  THAN  8,  SET  GA*8 

NOP 

G 

GAS 

STO 

GA 

G 

CLA 

GT 

G 

GA9 

SUB 

TEN 

ABSCISSA  DEC.  POINT  IS  MOD  10 

G 

TPL 

GA9 

G 

ADD 

TEN 

TZE 

GAIO 

G 

TPL 

GAIO 

G 

ZAC 

NEG.DEC.PT.  POSITION  «  0 

G 

GAIO 

STO 

GT 

G 

CLA 

SBV 

SBV,  COLUMNS  available  for  each  abscissa  VALUE 

STO 

G9 

G9  «  SBV  FIELD  WIDTH  FOR  ABSCISSA  VALUES 

G 

CAS 

GT 

IE  GT  GTR  THAN  OR  EUU  TO  G9,  SET  GT-G9-1 

TRA 

PASSA 

ENSURES  DEC  POINT  INSIDE  FIELD 

NOP 
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SUB  ONE 

STO 

G7 

C<J- 1 

G 

PASS* 

CLA  TWELV 

TWELVE  SPACES  ON  LEFT  FOR  ORDINATE 

AND 

LABEL 

AOO  G7 

STO  G6 

G6  =  07*12  FIELD  WIOTH,  LEFT 

ABSCISSA  VALUE 

CLA  REEL 

T 

SUB  TOTAL 

SUB  G6 

TPL 

PASS5 

reclt-total-G6.lt. 0, 

REOUCE 

G6 

G 

ADD  G6 

STO  G6 

G6  REOUCEO  TO  NUMBER  OF  COLUMNS 

AVAILABLE 

PASS5 

CLA  G6 

WHENEVER  G7.GE.G6,  G7  =  G6-1 

CAS 

G7 

ENSURES  DEC  POINT  INSIDE  LEFTMOST 

FIELD 

TRA 

EXIT 

NOP 

SUB 

ONE 

STO 

G7 

G6-  1 

G 

* 

G 

♦ 

G 

* 

SET  THE  FORMATS 

G 

* 

G 

EXIT 

LOU 

G3 

HJS 

TSL 

BCOCON 

HJS 

SLW 

FM  LA 

HJS 

CLA 

G3 

IF  G3  IS  NEGATIVE  SET  SCALE 

HJS 

TZE 

HJSl 

FACTOR  IN  FMl  TO  NEGATIVE 

HJS 

TPL 

HJSl 

HJS 

MSN 

FM  IA 

HJS 

HJSl 

LOO 

G4 

HJS 

TSL 

BCOCON 

G 

SLW 

FM  IB 

G 

LOO 

G5 

HJS 

TSL 

BCOCON 

HJS 

SLW 

FM3A 

HJS 

CLA 

G5 

IF  G5  IS  NEGATIVE  SET  SCALE 

HJS 

TZE 

HJS2 

FACTOR  IN  FM3  TO  NEGATIVE 

HJS 

TPL 

HJS2 

HJS 

MSN 

FM3A 

HJS 

HJS2 

LOO 

G6 

HJS 

TSL 

BCOCON 

G 

SLW 

FM3B 

G 

LDO 

G7 

G 

TSL 

BCOCON 

G 

SLW 

FM3C 

G 

SLW 

FM3F 

G 

LOO 

NV 

G 

TSL 

BCOCON 

G 

SLW 

FM3D 

G 

LOQ 

G9 

G 

TSL 

BCOCON 

G 

SLW 

FM3E 

G 

* 

G 

ZAC 

G 

RETURN 

PLOT  l 

EXIT,  SUCCESSFUL  PL0T1 

SPACE 

5 

G 
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* 


♦ 

PLOT2 

* 


* 


* 


* 


* 


* 


PLOT  2 

REM 

REM 

REM 

REM 


HJS 

MAIN  JOB  OF  PICIT2  IS  TO  REPEATEOLY  LAY  DOWN  SAMPLE 
GRIDLINE  IDASH  TO  DASH-WORDS*!)  FOLLOWED  BY  ISBH-l) 
NON-GRID  L  I '4E i  I  BLANK  TO  BLANK-WORDS* 1 1  TO  FORM  THE 
GRID  IN  THE  IMAGE  REGION 


SAVE  1,2 


ENTRY  TO  PLOT2 


ST  2 

WR0N3 

STZ 

WR0N2 

CLA 

WROfll 

TNZ 

GA12 

CLA 

3,4 

STA 

PLT22 

STA 

PLT23 

STA 

PLT37 

CLA*  4,4 

TSX 

TSTrP, 2 

TRA 

BAD 

STO 

XMAX 

CLA*  5,4 

TSX 

TSTFP.2 

TRA 

BAD 

STO 

XMIN 

CLA*  6,4 

TSX 

T5TFP, 2 

TRA 

BAD 

STO 

YMAX 

CLA*  7,4 

TSX 

TSTFP.2 

TRA 

BAD 

STO 

YMIN 

CIS 

XMIN 

FAD 

XMAX 

TZE 

BAD 

STO 

SPANX 

CLS 

YMIN 

FAD 

YMAX 

TIE 

BAD 

STO 

SPANY 

CLA 

LINES 

TSX 

FLOAT, 2 

FDP 

SPANY 

STO 

U 

CLA 

TOT 

TSX 

FLOAT, 2 

FOP 

SPANX 

STO 

V 

CLA 

NV 

WRON3  *  0  CLEAR  ERROR  FLAGt  MISSING  PLOT2 

CLEAR  ERROR  FLAG*  PLOT 2 


IMAGE  ADDRESS 


XMAX,  MAX.  ABSCISSA  VALUE 
TES1  FOR  FLOATING  POINT  ARGUMENT 


G 

G 


XMIN*  MIN.  ABSCISSA  VALUE 
TEST  FOR  FLOATING  POINT  ARGUMENT 


G 

G 


YMAX,  MAX.  ORDINATE  VALUE 
TEST  FOR  FLOATING  POINT  ARGUMENT 


G 

G 


YMIN*  MIN.  ORDINATE  VALUE 
TEST  FOR  FLOATING  POINT  ARGUMENT 


G 

G 


G 


ERROR  IF  XMIN  . 6Q-,  XMAX 
SPANX  *  XMAX  -  XMIN  ABSCISSA  SPAN 


G 

G 


ERROR  IF  YMIN  «EQ.  YMAX 
SPANY  •  YMAX  -  YMIN  ORDINATE  SPAN 


G 

G 


U  «  L INES/SPANY  NUMBER  OF  LINES  PER  UNIT  Y 

G 


V  *  TOT/SPANX  NUMBER  OF  COLUMNS  PER  UNIT  X 

G 


12  0  0  0  O  O  O  O  O 


-10S- 


TSX 

FLOAT, 2 

STO 

TEMP 

CIA 

SPANX 

FDP 

TEMP 

STQ 

DELTA 

DETLX  *  SPANX/NV  X  INCR  BETWN  VERT  GRID  LINES 

* 

CIA 

NH 

G 

TSX 

FLOAT, 2 

STO 

TEMP 

CLA 

SPANY 

FDP 

TEMP 

STQ 

DELTY 

DELTY  *  SPANY/NH  Y  INCR  BETWN  HORZ  GRID  LINES 

* 

STZ 

I 

1=0 

INITIALIZE  WORD  COUNTER  FOR  IMAGE  REGION 

G 

STZ 

J 

J=0 

INITIALIZE  LINE  COUNTER  FOR  IMAGE  REGION 

* 

G 

* 

G 

TNEXT 

STZ 

K 

K*0 

INITIALIZE  WORD  COUNTER  FOR  HORZ  GRID  LINE 

* 

G 

LNEXT 

CLA 

K 

LOOP  TO  PLACE  ONE  HORIZONTAL 

G 

PAX 

0,2 

GRID 

LINE  IMAGE  (DASH  REGION)  INTO  IMAGE  REGION 

ADD 

I 

PAC 

0,4 

G 

CLA 

DASH, 2 

FLT22 

STO 

♦  *.4 

IMAGE 

G 

CLA 

K 

ADD 

ONE 

STO 

K 

K  « 

K+ l  INCREMENT  WORD  COUNTER  FOR  LINE 

SUB 

WORDS 

TNZ 

LNEXT 

IF  NON-ZERO,  LINE  NOT  FINISHED 

G 

* 

G 

* 

CLA 

I 

G 

ADD 

WORDS 

STO 

I 

I  * 

I ♦WORDS  INCR  WORD  COUNTER  FOR  NEW  LINE 

* 

G 

CLA 

NH 

SEE  IF  FINISHED 

G 

SUB 

J 

TZE 

GA  13 

WHEN  J.EQ.NH,  IMAGE  GRID  COMPLETE 

G 

* 

CLA 

ONE 

G 

STO 

TEMP 

TEMP 

*  l  INITIALIZE  BETWEEN  GRID  LINE  COUNTER 

* 

G 

TFIN 

* 

STZ 

K 

7* 

II 

O 

INITIALIZE  WORD  COUNTER  FOR  EACH  LINE 

G 

LFIN 

CLA 

K 

LOOP  J  PLACE  SBH  NON-GRID  LINE 

G 

PAX 

0,2 

IMAGES  (BLANK  REGION)  INTO  THE  IMAGE  REGION 

ADD 

I 

PAC 

0,4 

G 

CLA 

BLANK, 2 

PLT23 

STO 

**  »  4 

IMAGE 

G 

CLA 

K 

ADD 

ONE 

STO 

K 

K  = 

X*1  INCREMENT  WORD  COUNTER  FOR  LINE 

SUB 

WORDS 
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TNZ 


LFIN 


* 

* 

* 

BAD 


* 

GA12 


* 

* 

GA13 


* 

* 


CLA  ! 

ADD  WORDS 
STO  I 
CLA  TEMP 
ADD  ONE 
STO  TEMF 
SUB  SBH 
TNZ  If  IN 

CLA  J 

ADO  ONE 
STO  J 

TRA  TNEXT 


IF  NON-ZERO,  LINE  NOT  FINISHED  G 

G 


I  *  I ♦WORDS  I NCR  WORD  COUNTER  FOR  NEW  LINE 

TEMP  *  TEMP* I  INCREMENT  BETWN  GRID  LINE  COUNTER 

IF  NON-ZERO.  WORE  LINES  REQUIRED  G 

G 


J  *  J*l 


INCREMENT  LINE  COUNT  FOR  IMAGE  REGION 
RETURN  FOR  ANOTHER  HORIZ.  GRID  LINE  G 

G 

G 

r 


RETURN  2,  UNSUCCESSFUL  PLQT2 
SE1  ERROR  FLAG ,  PLOT'! 


♦ 

PL0T3 

* 


CAL 

CTAPc 

CALL 

wcata 

PZE 

FORM 

PZE 

ER2.0, 1 

PZE 

WRONG, 0,3 

PZE 

0 

CLA 

ONE 

SET  *  NO  PLOT  2  *  FLAG 

STO 

WR0N3 

CLA 

FPTWU 

STO 

WRQN2 

RETURN 

PL0T2 

SPACE 

5 

PLOT  3 

REM 

PL0T3  EXAMINES 

THE  DATA  POINT  TO  MAKE  SURE  IT  I 

REM 

floating  po: 

NT 

AND  THcN  PLACES  I  IN  THE  PROPER 

REM 

IN  THL  IMAGE 

RP 

01  ON 

SAVE 

1,2 

PL0T3  ENTRY  POINT 

STZ  FLAG1  FLAGl  = 

0 

PLOT 3  RETURN  PRESET  TO  ZERO 

CLA 

3, A 

ADDRESS,  PLOTTING  CHARACTER 

STA  PL  T  36 

CLA 

A, 4 

BASE  ADDRESS,  X  COORDINATES 

STA  PIT 

35 

CLA 

5,4 

BASE  ADDRESS,  Y  COORDINATES 

XTA  PL  T  34 

CLA  WR0N1 

ORA  WRDN2 

TNZ 

GA  16 

OUT  IF  BAD  PLOT!  OR  PL0T2 

ORA 

WR0N3 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

HJS 


G 

G 


0 

G 

G 

G 


G 

G 

G 

G 


TZE 

GAIA 

OUT  IF  PREVIOUS  PL0T2 

G 

• 

G 

cal 

otape 

P1CT3  M/0  PLGT2 

G 

call 

WOAT  A 

G 

pze 

FORM 

G 

PZE 

ER3.0.3 

G 

PZE 

0 

G 

GA16 

CLA 

FTHRF 

G 

• 

G 

GA18 

RETURN 

PLOT  3 

G 

SPACE 

5 

G 

• 

G 

GAIA 

CLA» 

6, A 

NDATA.  NO.  OF  POINTS 

G 

TSX 

FIX 

.2 

TNZ 

GA1? 

G 

* 

G 

CL  S 

FTHRE 

IF  NOATA  =  0.  NO  DATA  POINTS-  RETURN  MINUS  THREE 

TRA 

GA18 

G 

• 

G 

GAIT 

STO 

NDATA 

NC.  OF  POINTS 

G 

STZ 

K 

A  =  U 

INITIALIZE  DATA  POINT  COUNTER 

* 

G 

ltend 

LAC 

K.l 

G 

* 

G 

PIT  34. 

CLS 

•*, 

l 

YOU 

Y  COORDINATE  OF  I X ♦ 1 ) TH  OAT*  POINT 

TSX 

TSTFP.Z 

TEST 

FOR  FLOATING  POINT  ARGUMENT 

TRA 

GA2 1 

G 

FAD 

Y  MAX 

LRS 

35 

FMP 

U 

T  PL 

GA  19 

G 

FS8 

NG5 

TRA 

GA20 

G 

GA  19 

F'O 

N05 

G 

GA20 

TSX 

FIX1.2 

G 

STO 

I 

I  =  ( YMAX-YlK) )*U  +08-  0-5,  LINE  INDEX  FOR  DATA 

PT 

TZE 

PL T  35 

Y  LIES  .ON  TOP  CRID  LINE 

G 

TMl 

GA2  1 

REJECT  Y  IF  ABOVE  TOP  GRID  LINE 

G 

SUB 

LINES 

TZE 

PLT35 

Y  LIES  ON  BOTTOM  GRID  LINE 

G 

*<• 

TPL 

GA21 

REJECT  Y  IF  BELOW  BOTTOM  GRID  LINE 

G 

* 

G 

PL  T 35 

CLA 

**, 

1 

X(K) 

X  COORDINATE  OF  (KM)TH  DATA  POINT 

TSX 

TS  TFP»  2 

TEST 

FOR  FLOATING  POINT  ARGUMENT 

TRA 

GA21 

G 

FSB 

XMIN 

LRS 

35 

FMP 

V 

TPL 

GA22 

G 

FSB 

NOS 

TRA 

GA2  3 

G 

GA  22 

FAD 

N05 

G 

GA2  3 

TSX 

F  I  X  1 . 2 

G 

STO 

J 

J=l X( K)-XMIN)*V  *0R-  0.5,  COL  INDEX  FOR  DATA 

PT 

TZE 

GA?A 

X  LIES  ON  LEFTMOST  GRID  LINE 

G 

TNI 

GA2 1 

REJECT  X  IF  LEFT  OF  GRID 

G 

SUB  TOT 

TZE 

GA  24 

X  LIES  ON  RIGHT  GRID  LINE 

G 

TPl 

GA21 

REJECT  X  IF  RIGHT  OF  GRID 

G 

GA24 

LOQ 

totls 

G 

MPY  1 

LLS 

35 

ADO  J 

STO  L 

L*TOTLS*I *J»  CHARACTER  POSITION  !N  INAGE 

REGION 

TSX 

PLACE, 4 

PLACE  BCD  IN  L-TH 

G 

PIT  36 

PZE 

*•  BCD 

G 

PZE 

L 

G 

PLT37 

PIE 

••  IMAGE 

G 

• 

G 

"MEMO 

CLA  K 

AOC  ONE 

STO  K 

R-K*  I 

INCREMENT  DATA  POINT  COUNTER 

SU8  NDATA 

TN2  LTENO  » 

IF  NONZERO.  MORE  DATA  POINTS  TO  BE  PLOTTEO 

« 

G 

CLA  FLA 

G1 

PLOT  3 

RETURN 

£ 

RETURN 

PLOT  3 

G 

GA2I 

CLS 

FTHRfc 

PLCT3  REJECTED  POINT 

U 

G 

STC 

FLAG1 

G 

THENU 

G 

SPACE  ' 

5 

G 

•PLOT* 

1  FPL0T4  1 

HJS 

REN 

PL0T4  OECONPOSES  THE  STRING  OF  CHARACTERS  IN  LABEL 

REN 

AND 

WRITES 

THE  CURRENT  GRAPH  ON  TAPE  OTAPE 

• 

G 

* 

G 

PLOT* 

SAVE 

1.2 

ENTRY  TO  PLCT4 

G 

• 

G 

• 

*034 

ASSURES  ALL  ARRAYS  ARE  STORED  FORWARD 

HJS 

FPL0T4 

SYN 

PLOT  4 

THEREFORE  BOTH  ENTRIES  ARE  THE 

SAME  G 

* 

G 

CLA  WXuKl 

ORA  WRON2 

TNZ 

GA26 

OUT  IF  BAD  PLOT l  OR  PLOT? 

G 

* 

G 

* 

G 

ORA 

WR0N3 

G 

TZE 

GA27 

OUT  IF  PREVIOUS  PL0T2 

G 

• 

G 

CAL 

OTAPE 

UNSUCCESSFUL  PL0T4 

G 

call 

MOAT  A 

G 

PZt 

FORP 

G 

PZE 

ER3.0.3 

NO  PREVIOUS  £L0T2 

G 

PIE 

0 

G 

* 

G 

GA26 

CLA 

FPFOR 

RETURN  4,  UNSUCCESSFUL  PL0T4 

G 

RETURN 

PL0T4 

G 

* 

G 

GA27 

CLA 

A, 4 

LABEL  BASE  ADDRESS 

G 
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PLT41 


Caz? 


* 

CHECK 


HJS3 


* 


AOOB 


STA 

PL  TAJ 

STA 

PL  T 42 

CLA 

PLT37 

INAGE  AODRESS 

G 

STA 

AODB 

G 

STA 

CA49 

G 

LXA 

WORDS, 2 

WORDS  PER  LINE 

G 

SXO 

AODB ,2 

G 

SXO 

GA49,2 

G 

CLA« 

3,4 

NCHAR,  NO.  OF  CHARACTERS  IN  LABEL 

G 

TSX 

FIX, 2 

ADD 

ONE 

PAX 

0,4 

NCHAR+l 

G 

AXT 

6,1 

SET  CHAR.  COUNT  FGR  LABEL  WORD 

HJS 

CLA 

YHAX 

SET  TOP  LINE 

G 

STO 

YAXIS 

ORDINATE  VALUE 

G 

r 

LDQ 

♦*  LABEL 

GET  FIRST  LABEL  WORD 

VP 

G 

STO 

LABEL 

CLA 

LINES 

STO 

FI  XV 

F I XV  *  LINcS  MAXIMUM  HORIZONTAL  LINE  INDEX 

CAL 

I  FONT 

ANA 

IF  OUR 

IF  IFOMT  *  4,5,6,  OR  7,  OELETE  BOTTOM  GRID  L 

INE 

TZE 

GA29 

r 

o 

CLA 

LINES 

SUB 

ONE 

STC 

FI  XV 

F I XV  =  L'INES-1,  MAX  LINE  INDEX  WITH  NO  BOTT 

LINE 

r; 

STZ 

I 

INITIALIZE  LINE  COUNT  FOR  IMAGE 

G 

AXT 

0,2 

INITIALIZE  WORD  COUNTER  FOR  IMAGE  REGION 

HJS 

G 

CLA 

F  I XV 

G 

SUB 

I 

TNI 

GA30  ' 

fIXV-I  NEGATIVE*  IMAGE  PRINT  COMPLETG 

TXH 

HJS3»4» 1 

IF  LABEL  HAS  BEEN  COMPLETELY 

VI 

HJS 

CAL 

BLNKK 

PRINTED,  OR  IF  NO  LABEL  IS  WANTED, 

HJS 

SLW 

LABEL 

SET  LABEL  TO  BLANK 

HJS 

ZAC 

HJS 

LOQ 

I 

DYP 

S8H 

TNZ 

SKIP 

IF  NON-INTEGRAL,  BYPASS  ORDINATE  PREPARATION 

CAL 

IFOMT 

ANA 

TWO 

TNZ 

SKIP 

IF  IFOMT*2,3,OR  6.  OELETE  OROINATE  VALUE 

CAL  OTAPE  GRID-LINE  IMAGE 

CALL  BOAT  A 

PZE  FMl 

PZE  LABEL « 0 t 1  LABEL  CHARACTER 

PZE  YAXIS.0,1 

PZE  **,<),*♦  IMAGE, 0, WORDS 

PZE  0 

CLA  YAXIS  ADJUST  ORDINATE  VALUE 


0000000000 
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fse 

delty 

G 

STO 

YAXIS 

G 

• 

G 

GA50 

CAL 

A008 

MOVE  TO  NEXT  IMAGE  LINE 

G 

A00 

WORDS 

G 

STA 

ADDB 

G 

STA 

GA49 

G 

CLA 

I 

COUNT  LINES 

G 

ADO 

ONE 

G 

STO 

I 

G 

• 

G 

TNX 

CHECK, 4,  l 

DECREMENT  LABEL  CHAR.  COUNT 

HJS 

CAL 

LABEL 

SET  NEXT  LABEL  CHARACTER 

G 

ALS 

6 

G 

TIX 

GA32 ,1*1  DECREMENT 

CHARACTER  COUNTER  IN  LABEL  WORD 

HJS 

AX  T 

6,1  REINITIALIZE  CHAR.  COUNT  FOR  LABEL  WORD 

HJS 

TXI 

♦  ♦ l »2  »-l 

MOVE  TO  AND 

G 

PLT42 

CAL 

**,2  LABEL 

GET  NEXT  LABEL  WORD 

G 

GA32 

SLW 

LABEL 

SAVE  LABEL  WORD 

G 

• 

G 

TRA 

CHECK 

AROUND  FOR  NEXT  LINE 

G 

• 

G 

* 

G 

SKIP 

CAL 

OTAPE 

WRITE  IMAGE,  NON-GRID-L iNE 

G 

CALL 

WOATA 

G 

PZE 

FM2 

G 

PZE 

LABEL, 0,1 

G 

GA49 

PZE 

♦*,0,**  IMAGE, 0, WORDS 

G 

PZE 

0 

G 

• 

G 

TRA 

GA50 

G 

• 

G 

* 

G 

GA30 

CAL 

IFOMT 

G 

ANA 

ONE 

TNZ 

EXIT2  IF  IFOMT* 

1,3,5,  OR  7,  DELETE  ABSCISSA  PRINTOUT 

* 

G 

CLA 

XHIN 

FIRST  ABSCISSA  VALUE 

G 

STO 

ABS 

G 

LXA 

NV ,  4 

FORM  ABSCISSA  VALUES 

G 

TXI 

**1,4,1 

NVM 

G 

SXO 

GA51 ,4 

r 

G 

TXI 

*♦1 ,4,-1 

NV 

G 

AXT 

0,1 

G 

CLA 

ABS,  l 

G 

LS2 

FAO 

DEL  T X 

G 

STO 

ABS* 1 , 1  *  ' 

G 

TXI 

*♦1, 1,-1 

G 

TIX 

L  S?  » 4 ,  l 

G 

• 

PUT  OUT  THE  ABSCISSA  LINE 

G 

CAL 

OTAPE 

CALL 

WDATA 

G 

PZt 

FM3 

G 

GAS  1 

P/r 

ABS,0-**  NV ♦ 1 

G 

PZE 

0 

G 
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6XIT2 


* 

* 

OMIT 


* 

GA37 


ZAC 

RETURN 

SPACE 

OMIT 

SAVE 

CLA* 


* 

* 

PLTAPE 


PLOT* 

5 


2 

3,* 


HJS 

HjS 

G 


TSX  TF I X.2 


DELETE  PRINTOUT  SEGMENTS 

ARGUMENT  TAKEN  MOD  8,  AS  FOLLOWS 
ARG=l,  OELETE  ABSCISSA  VALUE  PRINT 
ARG  *  2.  DELETE  ORDINATE  VALUE  PRINTOUT 


TNI 

GA37 

ARG=3» 

l  AND  2 

ORA 

I  FONT 

ARG**  t 

DELETE  BOTTOM  GRID  LINE 

SLW 

I  FONT 

RETURN 

OMIT 

ARG=5, 

l  AND  * 

COM 

ARG  =  6 1 

2  AND  * 

ANA 

IFOMT 

ARG=7, 

i»  2 »  AND  *» 

SLW 

IFCMT 

RETURN 

OMIT 

TO  RESTORE.  EAc  CUTE  CHIT 

REM 

SPACE 

PLTAPE 

SAVE 

CLA* 

TSX 

STO 

RETURN 

SPACE 

BCDCON 


2 

3.* 

FIX.2 

OTArE 

PLTAPE 

5 


WITH  THE  NEGATIVES  OF  THE  ABOVE  ARGUMENTS 


CHANGE  OUTPUT  TAPE  NO. 


CONVERT  INTEGER  IN  MQ 
TO  BCD  1999199  MAX.) 
IN  LOGICAL  AC. 


G 

G 

G 

G 

G 


BCDCON  AXT 

o 

* 

# 

« 

G 

AXT 

6»* 

G 

ZAC 

G 

BCD!  SLW 

8CD2 

G 

ZAC 

G 

DVP 

TEN 

G 

XEC 

BCD2  .  * 

G 

GRA 

BC02 

G 

T  I X 

BCDl. *» l 

G 

TRA* 

BCDCON 

G 

NOP 

U 

ALS 

6 

G 

ALS 

12 

G 

ALS 

18 

G 

ALS 

2* 

G 

ALS 

30 

G 

BC02  OEC 

0 

G 

* 

G 

000000  o  ooooooooo  oooo 
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SPACE 

5 

G 

• 

TSTFP 

TSTFP 

TZE 

2,2 

CHECK  FOR  FLOAT  1 NG  POINT 

G 

STO  TEMP 

ITS  CALLING  SEQUENCE  IS 

SSP 

TSX  TSTFP, 2 

SUB  MSOO* 

THE  ARGUMENT  IS  IN  THE  ACCUMULATOR 

TNI 

GAl  I 

G 

CLA  TEMP 

TRA 

2,2 

FLOATING-POINT  RETURN 

G 

GAll 

CLA 

TEMP 

G 

TRA 

1,2 

NON-FLOATING-POINT  RETURN 

G 

SPACE 

5 

G 

• 

• 

FLOAT 

FLOAT 

ORA  CONST 

FLOAT  FLOATS  A  NUMBER  KNOWN  TO  BE  AN  INTEGER 

FAD  CONST 

THE  CALLING  SEQUENCE  IS 

TRA  1,2 

TSX  FLOAT ,2  WITH  ARGUMENT  IN  ACCUMULATOR 

SPACE 

5 

G 

• 

• 

FIX 

* 

CONVERT  ARGUMENT  TO  FULL-WORD  INTEGEG 

FIX 

SSP 

POSITIVE  RESULT 

G 

TFIX 

LRS 

26 

SIGNED  RESULT 

G 

TNZ 

FIX2 

IF  NON-ZERO,  CONSIDERED  FLOATING 

G 

LLS 

26 

ALREADY  FIXED 

G 

TRA 

1,2 

FIX? 

LLS 

26 

RESTORE  FLOATING  NUMBER 

G 

F I X  l 

UFA  CONST 

FIXES  A  NUMBER  KNOWN  TO  BE  IN  FLOATING  POINT 

LRS  27 

ZAC 

G 

LLS  27 

TRA  1,2 

SPACE 

5 

G 

* 

• 

PLACF,  PLACB 

♦ 

PLACE  BCD  CHARACTER  IN 

G 

* 

I-TH  CHARACTER  POSITION  OF 

G 

• 

A  SPECIFIED  REGION 

G 

• 

TSX  PLACE, A 

G 

• 

PZE  BCD 

G 

• 

PZE  I 

G 

* 

PZE  REGION 

G 

PIACF 

MSM 

GA36 

-GA36,  FORWARD  ARRAYS 

G 

TRA 

GA35 

G 

• 

G 

PLACB 

MSP 

GA36 

♦  GA36 ,  BACKWARD  ARRAYS 

G 

* 

G 

GA35 

LOQ* 

2,  A 

I 

G 

STO 

TEM 

G 

ZAC 

G 

DVP  SIX 

STQ  TEKl 

TEHl*l/6,  INDEX  OF  WORD  CONTAINING  CHAR  POSITION 

LAC 

TEMU2 

SET  FOR  FORWARD  ARRAYS 

G 
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HIT 

GA36 

TEST  IF  FORWARD  ARRAY 

G 

LXA 

TEMl, 

2 

NO,  SET  FOR  BACKWARD  ARRAYS 

G 

HPY 

SIX 

G 

STO 

TEH1 

T£N1*TEM1*6,  CONTAINS  TOTAL  CHARACTER  POSITIONS 

CLA 

TEH 

TEH,  CONTAINS  THE  CHARACTER  POSITION 

SUB 

TEMl 

PAX 

0,1 

TEH- TEMl  =  CHARACTER  POSITION  IN  WORD 

CLA 

3,4 

STA 

GA4 

G 

STA 

GAS 

G 

CAL 

PLCHK,  l 

GA4 

ANA 

**,2 

ZERO  THE  CHARACTER  POSITION 

G 

SLK* 

GA4 

G 

CAL* 

1,4 

GET  AND 

G 

ANA 

MASK 

ISOLATE  THE  CHARACTER 

G 

TNX 

GA5, 1 

.0 

G 

GA6 

ARS 

6 

SHIFT  IT  INTO  POSITION 

G 

T I X 

GA6, 1 

,1 

G 

GA5 

ORA 

**,2 

PUT  CHARACTER  INTO  WORD 

G 

SLW* 

GA5 

G 

TRA 

4,4 

RETURN 

* 

G 

GA36 

PZE 

0 

-•FORWARD  ARRAY,  {-BACKWARD  ARRAY 

G 

SPACE  5 

G 

ABS 

BSS 

COL/ 2 

A8SCISSA  VALUES 

G 

BCI 

3, 

I 

BCI 

1.1 

SAMPLE  LINE  IMAGE 

G 

BCI 

1, 

FOR  NUN-GRID  LINES. 

BCI 

I,  I 

THIS  IS  LINE  IMAGE  USED  IN  STANDARD  GRID. 

BCI 

1, 

EXECUTION  OF  PLOT  1  SETS  UP  NEW  VALUES. 

BCI 

1. 

I 

BCI 

1,1 

BCI 

1, 

BCI 

l,  I 

BCI 

1, 

BCI 

!♦ 

1 

BCI 

1,1 

BCI 

1, 

BCI 

1,  I 

BCI 

1, 

BCI 

1, 

I 

BLANK 

BCI 

1,1 

FIRST  WORD  UE  LINE  IMAGE  FOR  NON-GRID  LINES 

8LNKK 

BCI 

1, 

CONST 

OCT 

233000000000 

* 

BCI 

3, 

- + 

BCI 

l»*~- 

— 

SAMPLE  LINE  IMAGE 

G 

BCI 

1, - 

— 

HORIZONTAL  GRID  LINES. 

BCI 

l,"  ♦ 

— 

THIS  IS  THE  LINE  IMAGE  USED  IN  STANDARD  GRID. 

BCI 

1, - 

— 

EXECUTION  OF  PLOTl  PRODUCES  NEW  VALUES. 

BCI 

It - 

BCI 

It* — 

— 

BCI 

if - 

— 

eci 

1,--* 

— 

BCI 

1.— 

— 

BCI 

1, - 

1 

i 

) 
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BCI 

It* - 

f  i 

BCI 

l. - 

BCI 

It— ♦— 

BCI 

1. - 

BCI 

1 

DASH 

BCI 

It* - 

1ST  WORD  OF  LINfc  IMAGE  FOR  HORIZ  GRID  LINES 

DEL  TX 

DEC 

0. 

INCR.  BETWEEN  VERTICAL  LINES 

G 

OELTY 

DEC 

0. 

INCR.  BETWEEN  HORIZ.  LINES 

G 

BC  1 

It- 

DSH  TO  DSH-5  USED  BY  PLOTl  TO  FILL  OUT  LAST 

BCI 

It— 

WORD  OF  HORIZ  LINE  IMAGE 

; 

BC  I 

I, — 

BCI 

l, - 

BCI 

it - 

OSH 

BCI 

It - 

EIGHT 

DEC  8 

ERl 

BCI 

ItOPLOTl 

r 

ER2 

BCI 

1 1 0PL0T2 

u 

ER3 

BCI 

3 1 ONO  PREVIOUS  PLQT2 

FIVE 

DEC 

5 

F I XV 

MAXIMUM  HORIZONTAL  LINE  INDEX  FOR  PRINTING 

flagi 

• 

* 

G 

♦ 

FORMATS  FOR  PRINTING  THE  IMAGE 

G 

fhi 

BCI 

It ( 1XA1, 

GRID  LINES 

G 

G 

FMl  A 

BCI 

It  0 

G3 

G 

BCI 

1 1  PF  9. 

r 

FMIB 

BCI 

1.  3 

G4 

Vj 

G 

# 

BCI 

2, , 1X20A6) 

G 

* 

NON-GRID  LINES 

G 

G 

FM2 

* 

BCI 

3 1 (  1XA1,  I0X20A6) 

G 

G 

♦ 

ABSCISSA  LINE 

G 

FH3 

BCI 

It ( 1 HO 

fl 

FM3A 

BCI 

It  0 

G5 

G 

BCI 

1 1  PF 

n, 

FM3B 

BCI 

it  15 

G6 

G 

BCI 

if* 

r 

FM3C 

BCI 

It  3 

G7 

G 

BCI 

Iff 

G 

FM10 

BCI 

It  10 

NV 

G 

BCI 

ltF 

G 

FH3E 

BCI 

It  10 

G9 

G 

BCI 

if* 

G 

FM3F 

BCI 

it  3 

G7 

G 

BCI 

it  ) 

G 

FORM 

BCI 

It (22A6) 

G 

FPONE 

DEC  I. 

FPTWO 

DEC  2. 

FTHRE 

DEC  3. 

* 

FPFOR 

DEC 

SPACE 

5 

C. 

* 

FORMAT  PARAMEfERS 

G 

* 

MODIFIABLE  BY  PLOT! 

G 
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*  0R01NATE  G3  P  f  9.G4 

*  LEFTMOST  ABSCISSA  G5  P  F  G6.G7 

*  OTHER  ABSCISSA  05  P  F  09.07 

G3  DEC  0 

G4  DEC  3 

05  DEC  0 

G7  DEC  3 

* 

G6  DEC  15  G7M2 

G9  DEC  10  SBV 

SPACE  5 

GW  I D  PZE  G  GRAPH  WIDTH  ( NUM8ER  OF  COLUMNS  FOR  IMAGE  ♦  1) 

I 

I E  YE  BCI  l.I 

I  PLUS  BCI  It* 

IFCMV  SWITCH  CONTAINING  INFORMATION  FROM  OMIT 

I  FOUR  DEC  4 

J 
K 
L 

LABEL 

LINES  DEC  50  ♦♦♦  MAXIMUM  HORIZONTAL  LINE  INDEX 

MS004  OCT  000777777777 
OCT  777777777700 
OCT  777777770077 
OCT  777777007777 
OCT  777700777777 
OCT  770077777777 


PLCMK 

OCT 

007777777777 

MASK 

OCT 

770000000000 

N05 

DEC 

0.5 

N999 

DEC 

.99 

NOATA 

'iUMBE  R 

OF 

DATA  POINTS  TO  BE  PLOTTED 

NH 

DEC 

5 

*** 

NUMBER 

OF 

HORIZONTAL  GRID  LINES 

NV 

DEC 

10 

*** 

NUMBER 

OF 

VERTICAL  GRiD  LINES 

ONE 

DEC 

l 

OTAPE 

DEC 

6 

OUTPUT  TAPE  ISET  BY  *  PLTAPE • ) 

RECLT 

PZE 

COL 

NUMBER 

CP 

COLUMNS  IN  OUTPUT  LiNE 

SBH 

DEC 

10 

*** 

NUMBER 

OF 

SPACES  BETWEEN  HORIZONTAL  GRIDLINES 

SBV 

PZE 

10 

*** 

SPACES 

BETWEEN  VERTICAL  GRID  LINES 

SIX 

DEC 

6 

5. 1 XF 

DEC 

6. 

SPAN* 

DEC 

0 

XMAX-XMIN 

SPANY 

DEC 

0 

ymax-ymin 

TEM 

TEMl 

TEMP 

TEN 

DEC 

10 

TOT 

DEC 

100 

*** 

MAXIMUM  COLUMN  INDEX 

TOTAL 

DEC 

101 

*** 

TOTAL 

COLUMNS  PER  LINE 

TOTL  S 

DEC 

102 

*** 

NUMBER 

OF 

BCD  CHARACTERS  PER  LINE 

TWELV 

DEC 

12 

TWO 

DEC 

2 

U 

DEC 

0. 

LINES  PER  UNIT  Y 

G 


G 


* 


G 

G 


oooooooooooo 
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V  DEC  0.  LINES  PER  UNIT  X  G 

WORDS  DEC  17  NUMBER  OF  MACHINE  LOCATIONS  PER  LINE 

WRONG  BCI  3*  IMPROPER  ARGUMENT 

WR0N1  OEC  0  *♦*  EQUALS  1  FOR  UNSUCCESSFUL  PLOU 

WR0N2  EQUALS  1  FOR  UNSUCCESSFUL  PL0T2 

WA0N3  OEC  l  ♦**  EQUALS  1  UNTIL  SUCCESSFUL  PL0T2 

YAXIS  G 

XMAX 

XNIN 

YMAX 

YMIN 

END  THIS  IS  THE  LAST  CARD 


$2  8MAP  10 

C  THIS  PROGRAM  IS  USED  IN  CONJUNCTION  WITH  THE  PLOTTING  ROUTINE 

C  TC  ALLOW  FORTRAN  INPUT-OUTPUT  FACILITY  IN  A  MACHINE  LANGUAGE 

C  PROGRAM 

«  - RANO  W038  I/O  ROUTINE  FOR  MAP  USERS 

*  WHO  WISH  TO  REFER  TO  THE  STD  FN4  I/O  PACKAGE, 

*  THIS  ROUTINE  IS  IOENTICAL  IN  FUNCTION 

*  TO  THE  7090  ROUTINE  X022.  THIS  7044  VERSION  IS 

*  OFFERED  COURTESY  OF  J  D  BABCOCK  WITH  THE  BLESSINGS 

*  C  W  ARMEROING  WHO  WAS  RESPOSIBLE  FOR  THIS  MESS 

*  ON  THE  7090. 

*  CALLS  ARE - 


CAL 

L  (LOGICAL 

TPAE  NO.  IN  DECREMENT) 

CALL 

(ROUTINE 

ENTRY) 

P/E 

FMT 

(BCD  ONLY) 

CP 

A1,T1,N1 

OP 

A2,T2»N2 

• 

•  •  • 

9 

•  •  • 

P/E 

0 

-  LAST  MUST  BE  ZERO. 

(RETURN) 

FMT  !S  LOCATION  OF  A  STANDARD  FN4  FORMAT  STATEMENT 
A I  *  T I  IS  THE  AODR.  OF  FIRST  DATA  WORO  (T»0,l,2> 

N I  IS  NUMBER  OF  WORDS  (AIDA-Ntli 
OP  IS  P/E  FOR  OIRECT  ADDRESSING 
OP  IS  M/E  FOR  INDIRECT  ADDRESSING 

ROUTINE  ENTRIES  ARE™ 

(1)  ROATA  -  BCD  INPUT 

(2)  IN  -  BCU  INPUT  FOR  STD  INPUT  UNiT 

(CAL  L  NOT  REQUIRED) 

(3)  WOATA  BCD  OUTPUT 

(4)  OUT  BCD  OUTPUT  FOR  STD  OUTPUT  UNIT 

ICAL  L  NOT  REQUIRED) 

(SI  PUNCH  BCD  OUTPUT  FOR  STD  PUNCH  UNIT 
(CAL  L  NOT  REQUIRED) 
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*  (61  WBIN  BINARY  OUTPUT 

*  (7)  RBIN  BINARY  INPUT 

* 


ENTRY 

MBIN 

ENTRY 

RBIN 

ENTRY 

IN 

ENTRY 

RDATA 

ENTRY 

OUT 

ENTRY 

MDATA 

ENTRY 

PUNCH 

* 

MB  I  N 

AXT 

0,0 

SXA 

IRA,  A 

AXT 

3,4 

TRA 

Dl 

♦ 

RBIN 

AXT 

0,C 

SXA 

IRA,* 

AXT 

2, A 

TRA 

Dl 

* 

PUNCH 

AXT 

0,0 

ALLON  CALLS  OF 

CAL 

*7 

TSL 

TRA 

MDATA 

OR 

* 

TSX  variety 

IN 

AXT 

0,0 

CAL 

*5 

STD  INPUT  UNIT 

RDATA 

AXT 

0,0 

SXA 

IRA, A 

AXT 

0  ,  A 

TRA 

Dl 

* 

OUT 

AXT 

0,0 

CAL 

*6 

STD  OUTPUT  UNIT 

MOATA 

AXT 

0,0 

SXA 

IRA, A 

AXT 

1 ,  A 

Dl 

SLW 

DT 

SAVE  LOGICAL  TAPE  NO, 

SXA 

ORA,* 

AXT 

7, A 

CAL* 

ANA 

ENTRY+7, A 

ADT 

TEST  FOR  TSX  OR  TSL  ENTRY 

TNZ 

07 

TlX 

*”3,4,1 

MSP 

ENTRY 

DR* 

AXT 

**  ,  A 

WAS  A  TRUE  TSX  ENTRY 

TRA 

D6 

07 

PAG 

0*  A 

FIX 

TXI 

*♦ i»4»— i 

UP  TRY  TO  SIMULATE 

CLA 

IRA 

TSX  ENTRY 

STA 

TSLRA 

SAVE  PROG.  IRA 

SXA 

IRA, A 

LXA 

ORA, A 

MSM 

ENTRY 
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06 

CL  A 

SEL  *4 

TAPE 

srn 

OSEL 

SET  UP  TSH,  STH 

CLA 

END, 4 

ANO  FIL,  RTN 

STO 

OFNO 

CLA 

CV 1  •  4 

STO 

0T1 

CLA 

DV2»4 

STO 

OT? 

TXH 

OB  IN* 4, 1 

TESf  IF  BCO  OR  BIN  CALL 

CLA 

MNL 

WAS  BCD 

STO 

CNVT 

SET  CONVERT  CALL 

LXA 

IRA, 4 

CAL 

1,4 

SET  UP  FORMAT 

TXI 

*♦1,4. -1 

BUMP  FOR  FIRST  OATA  CALL 

SXA 

IR4.4 

ANA 

ADT 

ORA 

BCDFM 

TRA 

GO 

Ob  KN 

CAL 

NOP 

FOR  BINARY  NOP  2,4 

LOQ 

BNL 

STO 

CNV1 

GO 

SLW 

01  NT 

CAL 

OT 

CALL  TG  SET  UP 

TSX 

UTVAR.,4 

FILE  NAME 

* 

GIVE  INITIAL  CALL 

OSEL 

*** 

♦  * 

TSX  TSHI0,4  OR  TSX  STHI0,4 

* 

OR  TSBIQ. , STBIO. 

OFILE 

PZE 

♦  * 

'ML  XX. 

DFMT 

PZE 

♦*,0,** 

FORMAT  FROM  CALL  IBCD-FORMAT, 

* 

NOW  SET  UP  BASIC 

LOOP — 

LXA 

I R4 ,4 

PREPARE  TO  PULL  OUT  ARGUMENTS 

02 

SXA 

IR4.4 

CLA 

1,4 

TZt 

DENO 

CHECK  IF  DONE 

POC 

0,4 

COUNT 

TXL 

IR4,4,0 

OUT  IF  NONE 

SXD 

DTST,4 

SET  TEST  FOR  NO.  Or-  WORDS 

TPL 

03 

CHECK  IF  INOIRECT 

ANA 

AOT 

YES,,  KILL  DECREMENT  ♦  PRFX. 

ORA 

DCLA 

SET  UP  CLA 

SLW 

*♦1 

*** 

** 

GETS  CLA  A, T  IF  INDIRECT 

03 

STA 

05 

SAVE  A 

ANA 

TAG 

PICK  UP 

ORA 

SXD4 

PROGRAMMER  INDEX  REGISTER 

SLW 

*♦1 

AND  PUT  THE  COMP.  OF  IT 

*** 

** 

LOC 

04,4 

INTO  IR4 

SXD 

04,4 

OS 

AXT 

**»4 

COMPUTE  ADDR  OF  A-T 

04 

TXI 

*♦1,4,4* 

SXA 

OT  1 ,4 

SET  ADDR. S  OF  PUTS 

SXA 

0T2,4 

ANO  GETS 

AXT 

0,4 

BEGIN  LOOP 
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on 

*** 

NOP  OR  CLA  DATA,  T 

CNVT 

♦  ** 

*♦ 

T*0(1)-N  (TSL  HNL I 0.  1 

DT2 

*** 

NOP  OR  STO  OAT  A , T 

* 

TXI 

♦♦1.4,-1 

NOP  =  AXT  0,0 

otst 

TXH 

DTI, 4,** 

-N  IN  DECR 

IR4 

AXT 

**,4 

a 

TXI 

02, 4,-1 

PEND 

*** 

0* 

FINAL  EXIT,  TSX  FILIC.,  4 

* 

OR  TSX  RTNIO.,4 

* 

AXT 

RESTORE  AXT 
7,4 

0,0  ENTRIES 

PXA 

0,0 

STA* 

ENTRY+7,4 

T I X 

*-1.4,1 

A 

LXA 

1R4,4 

w 

MIT 

ENTRY 

MAS  IT  TSX  OR  TSL 

TRA 

DENOi 

WAS  TSX  OK-EXIT  VIA  2,4 

TXI 

*♦1,4, -2 

WAS  TSL,  CALCULATE  RETURN 

PXA 

0,4 

PAC 

0,4 

ADDR 

SXA 

DENDl-l ,4 

LXA 

TSLR4.4 

RESTORE  PROG.  IR4 

TRA 

** 

RETURN  TO  MAIN  PROGRAM 

OENOl 

* 

TRA 

2,4 

♦ 

BCOFH 

MZE 

**,,FMTSC. 

BNL 

TSL 

BNL 10. 

HNL 

TSL 

HNL 10. 

OT 

PZE 

LOGICAL  TAPE 

apt 

OCT 

777777 

SAVL  AOOR  AND  TAG 

DCLA 

CLA 

**,0 

CLA  ORDER 

• 

* 

TABLE  OF  INITIAL  CALLS 

TSX 

STBIO.,4 

WBIN*IR4  OF  3 

TSX 

TSBIO.,4 

RB I N* I R4  OF  2 

TSX 

STHIO.,4 

WOATA  x  IR4  OF  l 

SEL 

A 

TSX 

TSHIO.,4 

RDATA  »  IR4  OF  0 

▼ 

TSX 

MLR  1 0. , 4 

WDUN 

TSX 

RLRIO.,4 

:<BUN 

TSX 

FILIO.,4 

WOATA 

END 

TSX 

RTNIO.,4 

ROATA 

* 

SET 

UP  HNL I 0.  , BNL 10.  LOOPS 

CLA 

** » 4 

3 

AXT 

0,0 

2 

CLA 

**»4 

1 

OVl 

* 

AXT 

0,0 

0 

AXT 

0,0 

3 

STO 

♦  *,4 

2 

AXT 

0,0 

1 

DV2 

STO 

♦*»4 

0 
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* 

♦ 


ENTRY 

PZE 

PUNCH 

PZE 

WOAf  A 

PZE 

OUT 

PZE 

ROATA 

PZE 

if: 

PZE 

WBIN 

PZE 

RB1N 

♦ 

TSLR4 

PZE 

♦♦ 

TAG 

OCT 

000000700000 

SXD4 

SXO 

04*0 

NOR 

EQU 

0V1 

* 

EXTERN 

UTVAR. 

EXTERN 

TSBIO. 

EXTERN 

STBIO. 

EXTERN 

BNLIO. 

EXTERN 

TSHIO. 

EXTERN 

STHIO. 

EXTERN 

HNL I0» 

EXTERN 

FMTSC. 

EXTERN 

RTNIO. 

EXTERN 

FILIO. 

EXTERN 

WLRIO. 

EXTERN 

END 

RLRIO. 

AIBMAP  UTV 

ENTRY 

UTVAR. 

EXTERN 

ERLOC. 

UTVAR.  SXA 

UTVX, 4 

SAVE  RETURN  INDEX 

LAC 

UTVX, 4 

SXA 

ERLOC., 4 

LXA 

UTVX, 4 

LAS 

NFIIES 

STOP  IF  LOGICAL  TAPE  NUMBER  EXCEEDS 

TRA 

USTOP 

NUMBER  OF  FILES  IN  TABLE. 

NOP 

PAC 

.4 

CLA 

IOU, 4 

PICKUP  ADDRESS  OF  FCB  POINTER 

PAX 

,  4 

TXL 

USTOP-2,4,0 

STOP  IF  UNIT  IS  UNDEFINED 

UTVX  AXT 

**,4 

RESTORE  RETURN  INDEX 

STO 

2,4 

SET  LOCATION  OF  FCB 

TRA 

l»4 

RETURN  TO  MAIN  PROGRAM 

LXA 

UTVX, 4 

CLA* 

“1,4 

RESTORE  UNIT  DESIGNATION 

USTOP  TSL 

FEXEM. 

ERROR,  ILLEGAL  UNIT  REQUESTED. 

PZE 

EXIT, ,3? 

NO  OPT IGNAL  RETURN 

♦INPUT-OUTPUT 

LOGICAL  UNIT  TABLE 

♦ADDITIONS  OR 

DELETIONS  SHOULD  BE 

MADE  BETWEEN  IOU  AND  NFILES 

IOU  PZE 

FiLOO. 
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NFILES 


PZE 

FIL01. 

PZE 

FIL02. 

PZE 

FIL03. 

PZE 

FIL04. 

PZE 

FIL05. 

PZE 

FIL06. 

PZE 

FIL07. 

PZE 

F I  LOB. 

****** 

PZE 

FIL09. 

****** 

PZE 

F IL 10. 

****** 

PZE 

FIL11. 

****** 

PZE 

F1L12. 

****** 

PZE 

FIL13. 

****** 

PZE 

F I L 14. 

****** 

PZE 

F 1 L 15. 

****** 

PZE 

♦-I0U-1 

EXTERN 

FILOO. 

****** 

EXTERN 

FIL01. 

****** 

EXTERN 

FIL02. 

****** 

EXTERN 

FIL03. 

****** 

EXTERN 

F1L04. 

****** 

EXTERN 

F I L05. 

****** 

EXTERN 

FIL06. 

****** 

EXVERN 

FIL07. 

****** 

EXTERN 

FIL08. 

****** 

EXTERN 

FIL09. 

****** 

EXTERN 

FIL 10. 

****** 

EXTERN 

FILll. 

****** 

EXTERN 

FIL12. 

****** 

EXTERN 

FIL13. 

****** 

EXTERN 

F I L 14. 

****** 

EXTERN 

FIL15. 

****** 

EXTERN 

;-EXEM.,EXIT 

END 
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